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13.  *as TRACT  (Mtjnmum  200  wom) 

)  In  the  course  of  performing  lethality  assessments  for  the  Defense  Nuclear  Agency  Single  Pulse  Laser  Lethality  and  Target 
Hardening  (LTH-3)  Program,  it  is  necessary  to  consistently  combine  many  sources  of  uncertainty  which  contribute  to  the 
overall  uncertainty  associated  with  algorithms  used  in  modelling  laser  interaction  and  target  response.  A  computer  code  called 
BETAFACT  was  developed  to  perform  this  task. 

A  Monte  Carlo  approach  is  used  in  BETAFACT  to  develop  an  algorithm  results  distribution  which  reflects  the  combined  effect 
of  all  sources  of  uncertainty  which  contribute  to  the  overall  algorithm  uncertainty.  Statistics  are  evaluated  for  the  algorithm 
results  distribution  to  quantify  the  overall  associated  uncertainty  and  to  identify  a  probability  distribution  which  models  the 
algorithm  results  distribution.  Normal,  lognormal.  Beta,  and  generalized  uniform  distributions  are  used  in  the  code  to  model 
uncertainties.  The  code  runs  interactively,  prompting  for  data  and  analysis  modifications  as  required.  It  is  necessary  for  the 
user  of  the  code  to  have  a  very  modest  proficiency  in  FORTRAN  programming  in  order  to  make  a  specific  algorithm  accessible 
to  the  code  via  user  provided  subroutines. 
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SUMMARY 


At  APTEK,  Inc.  we  use  an  algorithm  based  probabilistic  methodology  to  perform 
lethality  assessments  for  the  Defense  Nuclear  Agency  Single  Pulse  Laser  Lethality 
and  Target  Hardening  (LTH-3)  Program.  The  algorithms  utilized  relate  the  response 
of  laser  illuminated  targets  to  local  effects  (such  as  delivered  impulse)  induced  by  the 
impinging  laser  radiation.  In  general  the  algorithms  are  functions  of  laser  parameters 
(eg.,  fluence  and  pulse  duration)  and  target  design  details  (eg.,  physical  dimensions, 
construction  material  properties,  pressurization  levels,  etc.)  Typically  the  correlation 
between  algorithm  predictions  and  relevant  experimental  data  is  at  best  fair,  and  it  is 
always  the  case  that  target  design  details  are  imprecisely  known  since  we  are  dealing 
with  foreign  targets.  Thus  uncertainty  is  associated  with  our  algorithm  predictions 
because  of  the  uncertainty  associated  with  algorithm  parameters  compounded  with 
the  uncertainty  attributable  to  the  algorithms  themselves.  The  goal  of  the  work 
described  in  this  report  was  to  develop  and  implement  a  numerical  approach  for 
consistently  compounding  all  sources  of  uncertainty  in  order  to  obtain  a  quantified 
estimate  of  the  overall  uncertainty  which  should  be  associated  with  an  algorithm  and 
the  predictions  obtained  with  it.  The  combined  uncertainty  estimates  are  vital  inputs 
to  our  lethality  assessment  activities. 

When  one  is  attempting  to  consistently  combine  contributing  sources  of  error  or 
uncertainty  to  obtain  an  overall  uncertainty  estimate,  it  is  necessary  to  consider  sev¬ 
eral  important  factors.  Among  these  factors  are  the  characteristics  of  the  probability 
distribution  associated  with  each  source  of  uncertainty,  possible  correlation  between 
different  sources  of  uncertainty,  the  eventual  use  of  the  combined  uncertainty  esti¬ 
mate,  and  practical  aspects  of  implementing  a  procedure  to  obtain  the  combined 
uncertainty  estimates.  To  a  varying  extent  all  of  these  factors,  and  others  besides, 
influenced  the  development  of  the  finished  product  of  this  effort,  the  BETAFACT 
code. 

BETAFACT  is  a  Monte  Carlo  based  code  which  enables  the  user  to  numerically 
estimate  the  overall  uncertainty  associated  with  an  arbitrary  algorithm  which  has  an 
arbitrary  number  of  contributing  sources  of  uncertainty.  Normal,  lognormal,  Beta  and 
generalized  uniform  probability  distributions,  as  selected  by  the  analyst,  are  used  in 
the  code  to  model  contributing  uncertainties.  At  present,  for  a  given  analysis  all  un¬ 
certainty  sources  must  be  modeled  with  the  same  type  of  probability  distribution  (i.e., 
all  normal  or  all  Beta  distributed)  and  sources  of  uncertainty  must  be  uncorrelated. 
Since  the  code  runs  interactively,  once  an  algorithm  is  defined  it  is  an  easy  matter  to 
reanalyze  the  algorithm  for  different  probability  distributions  and  uncertainty  level 
specifications.  The  latter  capability  enables  the  user  to  perform  sensitivity  studies 
which  means  that  the  code  is  also  a  useful  decision  making  and  program  management 
tool. 

In  this  report  we  describe  the  technical  background  of  the  approach  used  to 
combine  uncertainty  sources,  discuss  implementation  of  the  theory  in  BETAFACT, 
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present  illustrative  examples  of  the  use  of  the  code,  and  provide  a  user’s  manual. 
Together  these  will  rapidly  acquaint  a  potential  user  with  the  operation  of  the  code. 
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CONVERSION  TABLE 


Conversion  factors  for  U.S.  Customary  to  metric  (SI)  units  of  measurement. 
MULTIPLY  — -  — *  BY  — ►  — ■  TO  GET 

TO  GET  4 —  4 —  BY  « —  < —  DIVIDE 


angstrom 

1.000  000  x  E  -10 

meters  (m) 

atmosphere  (normal) 

1.013  25  x  E  +2 

kilo  pascal  (kPa) 

bar 

1.000  000  x  E  +2 

kilo  pascal  (kPa) 

barn 

1.000  000  x  E  -28 

meter2  (m2) 

British  thermal  unit 

1.054  350  x  E  +3 

joule  (J) 

( thermochemical ) 

calorie  (thermochemical) 

4.184  000 

joule  (J) 

cal  (thermochemical)/cm2 

4.184  000  x  E-2 

mega  joule/m2  (MJ/m2) 

curie 

3.700  000  x  E  +1 

giga  becquerel  (GBq)* 

degree  (angle) 

1.745  329  x  E  -2 

radian  (rad) 

degree  Fahrenheit 

rK  =  (t°  f  +  459.67)/1.8 

degree  kelvin  (K) 

electron  volt 

1.602  19  x  E  -19 

joule  (J) 

erg 

1.000  000  x  E  -7 

joule  (J) 

erg/second 

1.000  000  x  E  -7 

watt  (W) 

foot 

3.048  000  x  E  -1 

meter  (m) 

foot-pound-force 

1.355  818 

joule  (J) 

gallon  (U.S.  liquid) 

3.785  412  x  E  -3 

meter3  (m3) 

inch 

2.540  000  x  E  -2 

meter  (m) 

jerk 

1.000  000  x  E  +9 

joule  (J) 

joule/kilogram  (J/kg) 

1.000  000 

Gray  (Gy)  ** 

(radiation  dose  absorbed) 

kilotons 

4.183 

terajoules 

kip  (1000  lu') 

1  448  222  x  E  t-3 

newton  (N) 

kip/inch2  (ksi) 

6.894  757  x  E  +3 

kilo  pascal  (kPa) 

ktap 

1.000  000  X  E  +2 

newton-  second/m2 
(N-s/m2) 

micron 

1.000  000  x  E-6 

meter  (m) 

mil 

2.540  000  x  E  -5 

meter  (m) 

mile  (international) 

1.609  344  x  F,  +3 

meter  (m) 

ounce 

2.834  952  x  E  -2 

kilogram  (kg) 

pound-force  (lbf  avoirdupois) 

4.448  222 

newton  (N) 

pound-force  inch 

1.129  848  x  E  -1 

newton-meter  (N*m) 

pound-force  inch 

1.751  268  x  E  +2 

newton/meter  (N/m) 

pound-force/foot2 

4.788  026  x  E  -2 

kilo  pascal  (kPa) 

pound-force/inch2  (psi) 

6.894  757 

kilo  pascal  (kPa) 

pound-mass  (lbm  avoirdupois) 

4.535  924  x  E  -1 

kilogram  (kg) 

pound-mass-foot2 

4.214  011  x  E  -2 

kilogram-  meter2 

(moment  of  inertia) 

(kg*m2) 

pound-mass/foot3 

1.601  846  x  E  +1 

kilogram/ meter3 
(kg/m3 

rad  (radiation  dose  absorbed) 

1.000  000  x  E  -2 

Gray  (Gy)** 

roentgen 

2.579  760  x  E  -4 

coulomb/kilogram 

(C/kg) 

shake 

1.000  000  x  E  -8 

second  (s) 

slug 

1.459  390  x  E  +1 

kilogram  (kg) 

torr  (mm  Hg,  0°C) 

1.333  22  x  E  -l 

kilo  pascal  (kPa) 

*  The  hecquerel  (Bq)  is  the  S!  unit  of  radioactivity;  1  Bq  =  1  event/s. 
**The  Gray  (Gy)  is  the  SI  unit  of  absorbed  radiation. 
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SECTION  1 
INTRODUCTION 


The  BETAFACT  code  was  developed  by  APTEK  to  support  the  preparation  of  yearly 
lethality  assessments  (APTEK,  1988  and  APTEK,  1990)  for  the  Defense  Nuclear 
Agency  (DNA)  Single  Pulse  Laser  Lethality  and  Target  Hardening  (LTH-3)  Program. 
The  results  reported  in  these  lethality  assessments  were  generated  using  an  APTEK 
modified  version  of  the  DNA  FAST  (Failure  Analysis  by  Statistical  Techniques)  code 
(Rowan,  1974).  Among  its  several  uses,  BETAFACT  is  utilized  to  generate  input 
required  in  the  modelling  and  execution  of  FAST  analyses. 

The  numerical  models  which  can  be  analyzed  with  FAST  are  defined  in  terms  of. 
mathematical  relations  called  algorithms.  A  typical  algorithm  consists  of  an  output 
quantity  which  is  a  function  of  several  input  parameters.  In  the  usual  case,  each  of 
the  input  parameters  does  not  have  a  precisely  known  value  but  does  have  a  most 
probable  (or  nominal)  value  and  a  distribution  of  possible  values  about  the  nominal. 
Such  a  parameter  is  termed  a  parameter  with  associated  uncertainty  in  this  report. 
It  is  usually  possible  (besides  being  convenient)  to  characterize  c  "'pproximate  the 
distribution  of  possible  values  for  a  parameter  with  associated  uncertainty  in  the  form 
of  standard  probability  distributions  such  as  the  normal  and  the  uniform  distributions. 

Clearly,  if  an  algorithm  is  a  function  of  one  or  more  parameters  with  associated 
uncertainty,  the  output  of  the  algorithm  must  also  have  associated  uncertainty  and  a 
distribution  of  possible  results.  The  output  uncertainty  and  output  distribution  are 
functions  of  the  input  parameter  uncertainties  and  distributions  and  the  details  of 
the  algorithm.  A  primary  function  of  BETAFACT  is  to  properly  combine  the  input 
parameter  uncertainties  and  distributions  in  order  to  determine  the  uncertainty  which 
should  be  associated  with  the  algorithm  output. 

An  algorithm,  viewed  from  an  overall  perspective,  may  also  have  an  associated 
uncertainty,  regardless  of  whether  or  not  its  input  parameters  have  associated  un¬ 
certainty.  This  situation  can  arise,  for  instance,  if  we  know  the  values  of  input 
parameters  accurately  (to  within  a  few  percent)  while  the  correlation  between  the 
parameters  provided  by  the  algorithm  is  only  known  approximately  (e.g.,  within  plus 
or  minus  50%).  In  this  situation,  the  input  parameter  uncertainties  could  be  viewed 
as  ignorable,  and  the  uncertainty  associated  with  the  algorithm  output  would  be  that 
applicable  to  the  algorithm  itself.  However,  if  the  input  parameter  uncertainties  are 
comparable  in  magnitude  to  the  algorithm  uncertainty  or  there  are  several  input  pa¬ 
rameters  with  smaller  but  still  not  ignorable  associated  uncertainty,  then  all  of  these 
sources  of  uncertainty  have  to  be  combined  to  obtain  the  uncertainty  which  should  be 
associated  with  the  algorithm  output.  The  BETAFACT  code  is  capable  of  handling 
both  of  these  situations. 

Since  the  code  combines  several  sources  of  uncertainty  into  an  overall  uncertainty 
associated  with  an  algorithm  output,  BETAFACT  can  also  be  used  to  perform  var¬ 
ious  types  of  sensitivity  studies.  For  instance,  suppose  one  has  an  algorithm  of  in- 
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terest  which  has  two  input  parameters,  each  with  significant  associated  uncertainty. 
By  varying  the  uncertainty  specified  with  each  of  the  two  parameters  and  observ¬ 
ing  the  effect  of  the  variations  on  the  overall  combined  uncertainty  estimated  with 
BETAFACT,  the  more  significant  of  the  two  uncertainty  sources  in  producing  the 
overall  uncertainty  can  be  identified.  While  such  a  procedure  may  not  be  necessary 
in  order  to  recognize  the  more  significant  uncertainty  source  in  a  simple  algorithm, 
recognition  of  that  source  can  be  much  more  difficult  in  a  complicated  algorithm.  As 
another  example,  again  suppose  an  algorithm  uses  two  parameters,  each  of  which  has 
associated  uncertainty  and  is  obtained  by  experiment.  A  project  manager  or  experi¬ 
menter  could  use  BETAFACT  to  evaluate  the  benefit  of  directing  program  resources 
to  reduction  of  the  parameter  uncertainties  both  in  terms  of  resources  expended  and 
reduction  in  overall  uncertainty  of  the  algorithm.  In  a  later  section  of  this  report  we 
will  provide  examples  illustrative  of  the  above  uses  of  BETAFACT. 

The  above  overview  provides  a  brief  description  of  the  function  of  BETAFACT 
and  some  simple  examples  of  how  the  program  can  be  used.  In  the  following  report 
section,  we  present  background  material  on  algorithms  and  treatment  of  their  sources 
of  uncertainty  and  review  some  properties  of  the  four  types  of  probability  distribu¬ 
tions  (normal,  lognormal,  Beta,  and  uniform)  which  are  available  in  BETAFACT  for 
modelling  parameter  and  algorithm  uncertainties. 

The  manner  in  which  algorithms,  probability  distributions,  and  uncertainties  are 
handled  numerically  in  BETAFACT  is  the  subject  of  Section  3.  This  section  is  par¬ 
ticularly  important  to  the  user  of  the  code  since  it  defines  the  K-factor  method  we 
employ  to  represent  uncertainty  specifications.  In  Section  4  we  consider  several  ex¬ 
ample  problems  which  illustrate  the  use  of  BETAFACT  and  verify,  in  part,  that  the 
theory  presented  in  Section  2  is  properly  implemented  in  the  code.  This  section  also 
addresses  the  method  which  must  be  used  to  get  a  user’s  algorithm  into  the  code. 

Section  5  of  this  report  is  the  User’s  Manual  for  BETAFACT.  Since  the  code  runs 
in  an  interactive  mode,  prompting  for  input  as  required,  the  user  will  likely  find  little 
need  for  the  User’s  Manual  after  the  program  has  once  been  used  successfully.  In 
Section  5  we  will  describe  the  actual  steps  required  to  setup  and  run  the  code  and 
explain  each  of  the  prompts  delivered  to  the  user  by  BETAFACT  during  execution. 
We  will  then  make  some  brief  concluding  remarks  in  Section  6.  Finally,  listings 
of  example  problem  algorithms  are  given  in  Appendix  A  and  a  complete  listing  of 
BETAFACT  is  presented  in  Appendix  B. 

Before  proceeding  to  Section  2,  we  make  note  that  BETAFACT  was  developed 
on  a  MicroVaxII  computer  operating  under  VMS.  The  code  is  written  in  standard 
FORTRAN-77  and  is  nearly  entirely  flexibly  dimensioned.  Although  not  verified,  it 
is  likely  the  program  will  also  execute  without  modification  on  personal  computers. 
Finally  we  mention  that  the  user  of  BETAFACT  must  have  simple  programming 
experience  with  FORTRAN  in  order  to  use  the  code.  This  is  a  requirement  because 
the  algorithm  of  interest  to  a  user  must  be  coded  in  a  FORTRAN  subroutine  so  that 
it  will  be  available  to  BETAFACT. 
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SECTION  2 

BACKGROUND  ON  ALGORITHMS  AND 
PROBABILITY  DISTRIBUTIONS 


The  primary  function  of  BETAFACT  is  to  determine  the  overall  uncertainty  asso¬ 
ciated  with  the  output  of  an  algorithm,  given  the  uncertainties  associated  with  the 
algorithm  and/or  its  input  parameters  and  corresponding  probability  distribution 
specifications.  The  code  performs  this  function  by  numerically  evaluating  the  algo¬ 
rithm  many  times  using  probabilistically  distributed  values  for  its  input  parameters, 
applying  (if  necessary)  algorithm  uncertainty  to  the  result  of  each  evaluation,  and 
computing  statistics  of  the  results  distribution  to  quantify  the  uncertainty  associated 
with  the  output.  This  can  be  done  for  an  algorithm  with  essentially  any  functional 
form  as  long  as  the  sources  of  uncertainty  associated  with  the  algorithm  can  be  mod¬ 
elled  adequately  using  either  normal,  lognormal,  Beta,  or  uniform  probability  distri¬ 
butions.  The  present  version  of  the  code  does  not  allow  different  sources  of  uncertainty 
in  an  algorithm  to  be  modelled  with  different  types  of  probability  distributions.  All 
uncertainties  associated  with  an  algorithm  must  be  modelled  as  normally  distributed 
or  Beta  distributed,  etc.  In  our  applications,  we  have  not  found  it  necessary  to  mix 
probability  distribution  types,  and  the  code  reflects  this  experience. 

The  purpose  of  the  first  subsection  (2.1)  of  this  report  section  is  to  provide  addi¬ 
tional  information  about  what  we  mean  by  the  concept  of  an  algorithm,  to  describe  a 
simple  method  for  quantifying  the  overall  uncertainty  associated  with  an  algorithm, 
and  to  give  an  explanation  of  why  numerical  treatment  of  algorithm  output  uncer¬ 
tainty  is  desirable  (and  in  some  cases  imperative).  In  subsection  2.2  we  then  present 
some  details  concerning  the  probability  distributions  used  in  BETAFACT  and  rele¬ 
vant  characteristics  of  the  distributions. 

2.1  ALGORITHMS  AND  UNCERTAINTIES. 

In  the  present  context,  an  algorithm  is  a  mathematical  function  which  relates  an 
output  or  response  (r)  to  one  or  more  input  parameters  (x,-,z  =  1,...,  number  of 
parameters  N).  Symbolically  we  write 

r  =  r(  xux2,...,xN)  (1) 

in  which  the  left  hand  side  is  the  output  of  the  algorithm  and  the  right  hand  side 
represents  the  functional  form  of  the  algorithm.  There  is  really  no  restriction  on  the 
form  of  the  function  on  the  right  hand  side  of  Equation  1  except  that  the  output  r 
must  not  be  an  argument  in  the  function  (i.e.,r  as  a  variable  must  not  appear  in  both 
sides  of  the  equation)  and  the  function  should  be  well-behaved  (i.e.  finitely  bounded) 
over  the  full  range  of  response  values  anticipated.  Otherwise,  the  function  r(x,-)  can 
be  polynomial,  exponential,  sinusoidal,  etc.,  and  involve  any  combinations  of  these 
functional  types. 
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A  specific  example  of  a  simple  algorithm  is  the  elementary  relation  between  the 
peak  bending  stress  a  in  a  straight  beam,  the  applied  bending  moment  M,  the  area 
moment  of  inertia  /  of  the  beam  cross  section,  and  the  extreme  fiber  distance  c  from 
the  beam  neutral  axis. 

Me 

<t=—  (2) 

In  this  form,  this  relation  is  an  algorithm  for  the  peak  bending  stress  <r.  Given  M,  c, 
and  /,  we  can  predict  the  peak  stress  in  the  beam.  However,  we  could  also  solve  this 
equation  for  M  in  terms  of  <r,  /,  and  c  to  obtain  an  algorithm  for  M.  If  we  accept 
the  assumption  that  the  formula  in  Equation  2  is  exact,  the  output  of  the  algorithm 
would  be  uncertain  only  if  some  of  the  input  parameters  had  associated  uncertainty. 

As  an  example  of  a  somewhat  more  complicated  algorithm,  we  consider  one  of¬ 
ten  used  in  the  LTH-3  Program;  the  so-called  HKL  algorithm.  One  version  of  this 
algorithm  is 

£  =  2  pah2  ^ 

It  is  used  to  predict  the  peak  axial  strain  e  which  results  in  a  metallic  cylinder  of 
thickness  h,  with  wall  material  density  p  and  ultimate  strength  cr,  which  is  loaded 
on  its  side  with  a  cosine  distributed  impulse  of  peak  intensity  I.  This  algorithm  is 
derived  based  on  theoretical  considerations.  Given  good  (i.e.,  known  or  measured 
with  low  associated  uncertainty)  values  for  the  parameters  on  the  right  hand  side 
of  the  equation,  the  strain  predictions  obtained  with  Equation  3  have  been  shown 
to  agree  to  within  about  ±50%  with  experimentally  measured  strains.  It  is  thus 
an  example  of  an  algorithm  which  has  an  associated  algorithmic  uncertainty.  If  we 
used  this  algorithm  to  predict  the  axial  strain  resulting  in  a  cylinder  whose  p  and  a 
values  were  inaccurately  known  but  we  knew  h  and  I  with  considerable  precision,  we 
would  have  the  situation  of  an  algorithm  with  uncertainty  associated  with  two  input 
parameters  and  the  algorithm  itself.  All  three  sources  of  uncertainty  contribute  to 
the  uncertainty  associated  with  the  strain  estimate  for  that  cylinder. 

Given  an  algorithm  of  the  form  of  Equation  1,  a  simple  standard  method  is  avail¬ 
able  for  estimating  the  uncertainty  associated  with  the  algorithm  output  if  one  or 
more  of  its  input  parameters  has  associated  uncertainty.  First,  let  Ui  represent  the 
uncertainty  associated  with  the  i-th  parameter  (x*)  in  the  algorithm.  We  assume 
that  U,  is  not  correlated  with  the  uncertainty  Uj  associated  with  the  j-th  parameter 
for  all  pairs  i  and  j.  (In  this  report  and  in  the  BETAFACT  code,  the  uncertain¬ 
ties  associated  with  algorithm  input  parameters  are  all  assumed  to  be  uncorrelated.) 
Typically  Ui  is  selected  or  specified  such  that  95%  of  the  possible  values  of  x,  fall  in 
the  range  associated  with  U ,■  about  the  nominal  (or  50%  or  best  estimate)  value  for  x,. 
(Relating  the  nominal  value  of  x,,  Ui  and  the  applicable  range  for  x,  are  considered 
in  greater  detail  in  subsection  2.2.)  Now  if  only  the  i-th  parameter  of  the  algorithm 
is  assumed  to  have  associated  uncertainty,  then  the  corresponding  uncertainty  Ur,  in 
the  algorithm  output  can  be  evaluated  as 

Ur,  =  ~Ui  (4) 

OXi 
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The  usual  situation  is  that  several  input  parameters  in  the  algorithm  have  associ¬ 
ated  uncertainty.  Equation  4  can  be  evaluated  for  each  of  these  uncertainty  sources 
considered  separately,  but  the  question  then  becomes  one  of  how  to  combine  the 
individual  contributions  to  obtain  a  reasonable  estimate  for  the  algorithm  output 
uncertainty.  A  recommended  combination  rule  (Coleman  and  Steele,  1989)  is  to  use 
the  root-sum-square  (RSS)  of  the  C/r,.  Thus,  if  we  represent  the  combined  output 
uncertainty  as  Ur,  the  RSS  expression  for  Ur  in  terms  of  the  UTi  is 

N  «  fa.  2 

^  =  =  (5) 

i=i  1=1  LL,X*  . 

Given  the  algorithm  r  and  the  parameter  uncertainty  specifications  £/;,  Equation  5  is 
easily  evaluated  to  yield  UT. 

It  is  sometimes  found  that  a  more  convenient  form  of  Equation  5  is  to  divide  both 
sides  of  it  by  r2  (i.e.,  the  algorithm  squared)  to  obtain 


This  equation  often  will  give  algebraically  simpler  expressions  for  U2  than  will  Equa¬ 
tion  5. 

Although  the  above  results  are  straightforward,  it  is  perhaps  useful  to  apply  them 
to  consideration  of  an  example  algorithm;  a  hypothetical  one  in  this  case.  Suppose 
r  =  A2B/C  for  which  the  nominal  parameter  values  are  A  =  10,  B  =  20,  and  C  = 
50  and  the  associated  uncertainties  are  Ua  =  ±2,  Ub  =  ±4,  and  Uc  =  ±15.  The 
nominal  output  of  the  algorithm  is  40.  Differentiating  the  algorithm,  we  find  that 


~—U 
r  dA  A 


l°LUg=V»  ~—Uc  =  — 

rdB B  B  rdC  C  C 


(7) 


and  thus 


Substituting  for  the  given  values  we  find  Ur  =  ±21.5.  If  we  set  Ug  =  Uc  =  0,  then  we 
find  that  UT  =  ±16  which  shows  that  Ua  accounts  for  the  most  significant  portion  of 
the  combined  uncertainty  even  though  the  uncertainty  associated  with  A(±20%)  is 
the  same  as  that  associated  with  B,  but  less  than  that  associated  with  C(±30%).  The 
exponent  of  A  in  the  algorithm,  which  has  an  absolute  value  greater  than  1,  is  seen 
to  enhance  the  contribution  of  Ua  to  the  overall  uncertainty.  Of  course  the  opposite 
is  also  true;  if  the  absolute  value  of  the  exponent  were  less  than  1,  the  contribution 
of  U a  would  be  lessened. 

The  above  method  for  estimating  UT  does  not  explicitly  account  for  the  actual 
probability  distributions  of  the  possible  values  of  the  parameters  with  associated  un¬ 
certainty.  Thus  the  UT  estimates  obtained  in  this  fashion  are  often  faulty.  In  the 
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special  case  that  these  parameter  values  are  normally  distributed,  it  can  be  shown 
rigorously  that  the  above  expression  gives  the  correct  value  for  UT ,  at  least  for  algo¬ 
rithms  of  the  form  r  —  Kx\x\x%....  in  which  if  is  a  constant  and  a,  b,  c, ...  are  also 
constants  (positive  or  negative).  For  other  probability  distributions  for  the  possible 
parameter  values  and  for  more  complicated  algorithms,  the  Ur  estimates  obtained 
from  Equations  5  or  6  may  or  may  not  compare  well  with  the  true  UT-  Numerical 
evaluation  of  an  algorithm  using  properly  distributed  values  of  the  uncertain  input 
parameters  and  statistical  analysis  of  the  distribution  of  evaluated  algorithm  results 
together  appear  to  be  the  only  method  available  for  obtaining  good  estimates  for 
the  combined  uncertainty  UT  associated  with  the  algorithm  output.  This  is  the  ap¬ 
proach  used  in  BETAFACT.  The  actual  implementation  of  the  approach  is  described 
in  Section  3. 

Numerical  evaluation  of  an  algorithm  allows  one  to  model  its  associated  uncer¬ 
tainty  sources  using  essentially  any  distribution  type  desired.  We  have  found  it  useful 
to  model  parameter  and  algorithm  uncertainties  with  normal,  lognormal,  Beta,  and 
uniform  probability  distributions.  The  forms  and  specific  characteristics  of  these 
distribution  types  are  summarized  in  the  following  subsection. 

2.2  PROBABILITY  DISTRIBUTIONS. 

The  distribution  of  possible  values  of  each  parameter  or  algorithm  with  associated 
uncertainty  is  modelled  with  a  normal,  lognormal,  Beta,  or  uniform  probability  dis¬ 
tribution  in  BETAFACT.  Each  of  these  probability  distribution  types  has  a  more  or 
less  convenient  mathematical  representation.  Before  presenting  these  functions  and 
summarizing  some  of  their  important  characteristics,  we  very  briefly  review  for  the 
reader  the  meaning  of  a  probability  distribution  and  describe  several  measures  gener¬ 
ally  used  to  characterize  a  probability  distribution.  We  will  see  in  Section  3  how  the 
latter  are  related  to  a  given  specification  of  uncertainty  for  a  parameter  or  algorithm. 

Let  f(x)  represent  the  probability  distribution  associated  with  the  distributed 
variable  x.  Then  f(x)dx  is  the  probability  that  x  has  a  value  in  the  interval  between 
x  and  x  +  dx.  Suppose  we  know  variable  x  only  takes  on  values  in  the  range  a  <  x  <  b. 
Since  x  must  fall  in  this  range,  the  integral  of  /(x)  over  the  full  range  of  possible  values 
is  unity,  provided  /(x)  is  normalized  (which  is  the  case  for  all  four  distributions  we 
will  consider  below).  Thus 


/  f(x)dx  =  1  (9) 

Ja 

For  some  distributions,  such  as  the  uniform  distribution,  the  limits  a  and  b  on  the 
range  of  x  variable  values  are  bounded;  for  others,  such  as  the  normal  distribution, 
they  are  not. 

The  mean  p  (or  expectation  value  E(x))  of  a  probability  distribution  is  defined 
by  the  integral 

p  =  E(x)  =  f  xf(x)dx  (10) 

J  a 
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(11) 


and  the  variance,  E(x  —  /z)2,  by  the  integral 

E(x  —  /z)2  =  f  (x  —  fi)2f(x)dx 
Ja 

From  the  variance  we  obtain  the  standard  deviation  cr 

a  =  y/E( x  -  n)2  (12) 

Now  suppose  we  have  collected  a  set  of  N  values  for  the  variable  x.  For  instance, 
if  x  is  a  physical  parameter,  the  set  of  values  could  be  the  measured  results  for  x 
as  found  in  many  experiments.  In  the  BETAFACT  context,  the  set  of  values  of 
interest  is  the  set  of  results  obtained  in  the  evaluation  of  a  algorithm  many  times. 
Let  us  further  suppose  that  the  set  of  x  values  have  been  sorted  in  ascending  order 
(smallest  to  largest),  the  range  from  the  smallest  to  largest  value  has  been  divided 
into  M  subintervals,  and  the  number  of  times  a  value  of  the  variable  occurs  in  each 
subinterval  has  been  determined.  Let  x,-  be  the  variable  value  corresponding  to  the 
midpoint  of  the  z-th  subinterval  and  /,  be  the  number  (or  frequency  of  occurrence) 
of  x  values  in  the  z-th  subinterval.  Clearly  N  =  Y/iLi  ft-  Also  if  we  plot  /,  as  a 
function  of  x,,  we  obtain  a  graphical  representation  of  the  probability  distribution  for 
the  variable  x.  Such  a  representation  is  often  termed  a  histogram.  If  we  have  collected 
a  sufficient  number  of  data  points  and  have  used  a  reasonable  number  of  subintervals 
in  segregating  the  data  into  frequency  of  occurrence  bins,  then  we  might  recognize 
the  plotted  distribution  as  being  representable  by  one  of  the  standard  mathematical 
probability  distributions. 

The  availability  of  the  M  pairs  (x„  /,)  enables  us  to  determine  a  mathematical 
distribution  which  more  or  less  closely  replicates  the  actual  distribution  of  variable 
values  once  we  have  selected  a  type  of  probability  distribution  to  model  the  variable 
distribution.  First  we  estimate  the  mean  x  of  the  variable  distribution 

i  A i 

*  =  Tr£*«7i  (13) 

iv  1=1 

and  the  variance  s2 

i  M 

s2  =  T7  £(xi  -  x)2fi  (14) 

JV  1=1 

The  evaluated  values  x  and  s  are  then  used  to  fit  the  appropriate  mathematical 
probability  distribution  to  the  distribution  of  x  variable  results  using  the  relations 
between  the  mean  and  standard  deviation  and  the  corresponding  distribution  given 
in  the  following  subsections  for  each  distribution  type  in  BETAFACT.  We  now  turn 
to  a  brief  description  of  the  general  form  and  characteristics  of  these  distributions. 


2.2.1  Normal  Distribution. 

The  functional  form  of  the  normal  or  Gaussian  probability  distribution  is 

(x~H)2' 


fix)  = 


1 


ay/ 


exp 


2<r2 


(15) 
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Figure  1.  Example  normal  probability  distributions  with  constant  mean  and  vari¬ 
able  standard  deviation. 


The  mean  of  this  distribution  is  p,  the  variance  is  a2,  and  the  standard  deviation 
is  cr.  The  range  of  possible  variable  values  is  unbounded  (— oo  <  x  <  oo).  The 
distribution  is  centered  on  and  symmetric  about  x  —  /z.  If  a  variable  x  is  normally 
distributed  with  mean  /z  and  standard  deviation  cr,  a  randomly  selected  value  for  x 
will  fall  95%  of  the  time  in  the  range  ( /z  —  1.96cr  <  x  <  p  4-  1.96<r).  Only  2.5% 
of  the  time  will  the  randomly  selected  value  for  x  be  below  this  range  and  2.5%  of 
the  time  above  this  range.  Finally  if  we  compute  x(~  /z)  from  Equation  13  and 
x(~  a)  from  Equation  14  for  a  set  of  variable  values  which  appear  to  be  normally 
distributed  and  then  substitute  these  values  into  Equation  15,  we  obtain  the  normal 
distribution  which  will  approximately  fit  the  distribution  of  variable  values.  As  an 
aid  in  visualizing  the  normal  distribution  and  the  effect  on  the  distribution  as  a  is 
varied  (with  the  mean  fixed),  Figure  1  shows  3  normal  distributions  for  3  different 
values  of  a  and  a  constant  value  for  /z. 

2.2.2  Lognormal  Distribution. 

Suppose  variable  x  has  a  probability  distribution  of  possible  values  which  has  a  mean 
M.  Further  suppose  the  distribution  for  z  is  characterized  by  the  existence  of  a 
positive  constant  K(>  1.0)  which  has  the  property  such  that  95%  of  the  time  a 
randomly  selected  value  for  the  variable  z  falls  in  the  range  M / K  <  r  £  K  x  M 
and  equally  as  often  above  M  as  below  M.  Then  the  distribution  of  possible  values 
for  z  is  lognormally  distributed  if  the  distribution  obtained  by  the  change  of  variable 
y  =  lnz  is  a  normal  distribution  in  y  with  mean  /z  =  In  M  and  standard  deviation 
a  =  lnAyi-96.  In  the  BETAFACT  (and  FAST)  context,  K  is  recognized  as  the 
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Figure  2.  Example  lognormal  probability  distributions  with  constant  mean  and  vari¬ 
able  standard  deviation. 

K-factor  which  characterizes  (or  quantifies)  the  uncertainty  associated  with  variable 
x. 


The  typical  form  of  the  lognormal  distribution  is  exhibited  by  the  three  example 
distributions  shown  in  Figure  2.  We  see  from  this  figure  that  the  effect  of  increasing  a 
or  K  on  the  lognormal  distribution  is  to  shift  the  distribution  peak  to  lower  variable 
values  and  significantly  stretch  out  the  high  end  tail  of  the  distribution. 


2.2.3  Beta  Distribution. 


The  Beta  probability  distribution  is  defined  on  the  unit  interval  0  <  x  <  1  by  the 
two  parameter  function 


/(*) 


xa~\l  -  x)0~l 
B(a,0) 


(16) 


in  which  the  normalization  factor  B(a,  0)  is  given  by 


B{a,0)  =  [' xa~'{l -xf~xdx  (17) 

Jo 

In  Equation  16,  a  and  0  are  two  parameters  whose  values  influence  the  shape  of  the 
corresponding  probability  distribution.  Both  a  and  0  are  restricted  to  be  greater 
than  zero,  and  in  fact  in  our  experience  with  BETAFACT,  all  the  Beta  distributions 
we  have  used  have  had  both  a  and  0  great  than  1. 

Figure  3  shows  the  Beta  distributions  obtained  using  3  sets  of  a,  0  pairs.  From 
the  figure  we  see  that  if  a  =  0,  the  resulting  probability  distribution  is  symmetric  and 
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Figure  3.  Example  Beta  probability  distributions  for  different  (a,  0)  parameter  val¬ 
ues. 

looks  somewhat  like  a  normal  distribution.  If  a  <  0  (and  a  >  1.0),  the  distribution  is 
skewed  to  the  left  and  the  opposite  is  true  if  a>  0  (and  0  >  1.0).  In  the  former  case, 
the  distribution  has  an  appearance  similar  to  a  lognormal  distribution.  The  ability 
to  skew  the  Beta  probability  distribution  to  the  right  by  proper  choice  of  a  and  0  is 
a  feature  of  the  function  not  shared  by  either  the  normal  or  lognormal  distribution. 

The  mean  m(=  y)  of  the  Beta  distribution  with  parameter  values  a  and  0  is 


m  = 


a 

a  4-  0 


and  the  variance  is 


a/3 

(a  +  0)2(a  +  0  +  1) 


Equations  18  and  19  can  be  inverted  to  give 


a  =  m 


0  =  (1  -  m) 


(18) 

(19) 


(20) 


However,  since  the  standard  deviation  of  the  distribution  must  be  positive,  the  fol¬ 
lowing  condition  must  also  be  satisfied  (Simons,  D.A.,  1988): 


a  <  min 


m 


(21) 


The  implications  of  the  above  restriction  on  a  will  be  considered  in  the  section  of  this 
report  in  which  implementation  of  Beta  distributions  in  BETAFACT  is  discussed. 
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If  we  have  a  set  of  variable  results  which  we  think  can  be  represented  with  a  Beta 
distribution,  Equations  16  and  20  indicate  how  this  might  be  done.  First  we  evaluate 
the  mean  m  and  the  variance  a2  of  the  set  of  variable  results  and  use  Equation  20  to 
evaluate  corresponding  values  for  a  and  (3.  Then  we  evaluate  Equation  16  to  obtain 
the  desired  distribution.  This  is  one  approach  used  in  BETAFACT.  A  somewhat  less 
demanding  (and  less  accurate)  method  is  also  implemented  in  the  code.  This  method 
will  be  described  in  Section  3. 


2.2.4  Uniform  Distribution. 

A  generalized  form  of  the  uniform  probability  distribution  is  used  in 
BETAFACT.  The  function  describing  this  distribution  is  defined  on  the  interval 
a  <  x  <  b  as 

/(*)  =  { ?  v:  “ -  *  <7  (22) 

l  2  b-m  m  —  X  —  " 

The  parameter  m  corresponds  to  the  value  of  x  such  that  50%  of  the  time  the  value 
of  x  is  less  than  m.  The  mean  \i  of  the  distribution  is 

a  +  b  +  2m 


and  the  variance  a2  is 


<x2  —  ~  [a2  +  62  +  m(a  +  b)  +  2m2]  —  j-  [a  +  b  +  2m]2 


For  m  =  i(a  +  b)  these  reduce  to 


P  =  \(a  +  b)  ^  ^  (b  ~  a  f 


as  required  for  a  symmetric  uniform  distribution.  Figure  4  shows  examples  of  the 
generalized  uniform  probability  distribution  associated  with  three  different  values  of 
m  for  fixed  a  and  b.  Clearly,  given  a  set  of  variable  values  which  can  be  modelled 
with  a  generalized  uniform  distribution,  we  need  to  determine  only  a,  6,  and  m  for 
the  set  of  variables  to  fit  the  required  uniform  distribution. 
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Figure  4.  Example  generalized  uniform  probability  distributions. 
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SECTION  3 
IMPLEMENTATION 


Four  types  of  probability  distributions  (normal,  lognormal,  Beta,  and  uniform)  are 
available  in  BETAFACT  for  modelling  the  distributions  of  parameters  and  algorithms 
which  have  associated  uncertainty.  Some  important  characteristics  of  each  of  these 
distributions  are  summarized  in  Section  2  of  this  report.  The  purpose  of  this  sec¬ 
tion  is  to  describe  the  specifics  of  how  the  code  uses  the  probability  distributions  to 
model  uncertainties,  and  the  mechanics  involved  in  developing  the  algorithm  results 
distribution  which  ultimately  provides  the  information  needed  to  quantify  the  overall 
uncertainty  associated  with  an  algorithm  output. 

Several  steps  are  involved  in  the  processing  of  an  algorithm  with  BETAFACT. 
These  steps,  with  minor  modifications,  need  to  be  performed  for  any  algorithm  an¬ 
alyzed  with  the  code.  The  steps  are  briefly  described  in  subsection  3.1.  The  precise 
manner  in  which  some  of  the  steps  are  implemented  in  BETAFACT  is  dependent 
on  the  probability  distribution  type  used  to  model  the  uncertainties.  Details  of  the 
implementation  are  described,  in  turn,  for  each  of  the  4  probability  distribution  types 
in  subsections  3.2  through  3.5. 

3.1  ANALYSIS  PROCESS  OVERVIEW. 

In  the  following  we  assume  that  the  algorithm  to  be  analyzed  with  BETAFACT  is 
reduced  to  the  form  of  Equation  1;  i.e.  the  algorithm  output  is  expressed  strictly  in 
terms  of  input  parameters,  any  number  of  which  have  associated  uncertainty.  It  has 
been  our  experience  that  at  least  one  parameter  in  the  algorithm  may  be  regarded  as 
an  independent  parameter  which  does  not  have  associated  uncertainty.  BETAFACT 
expects  there  to  be  at  least  one  such  parameter.  If  there  isn't  explicitly  such  a 
parameter,  one  can  easily  be  introduced  as  a  multiplier,  with  value  of  unity,  which 
acts  on  the  entire  right  hand  side  of  Equation  1. 

To  be  available  to  BETAFACT,  the  algorithm  must  be  coded  in  FORTRAN  in  a 
subroutine  (named  TF)  called  by  the  code.  This  subroutine  must  be  compiled  and 
linked  with  the  BETAFACT  object  file  to  obtain  the  algorithm  specific  executable 
form  of  the  program.  Details  of  how  this  is  accomplished  are  given  in  Section  5  of 
this  report. 

We  have  found  on  occasion  that  some  parameters,  either  with  or  without  uncer¬ 
tainty,  which  are  used  in  an  algorithm  are  most  conveniently  obtained  from  subordi¬ 
nate  algorithms.  For  instance,  the  thickness  parameter  h  in  the  algorithm  of  Equation 
3  may  be  a  known  function  of  the  impulse  intensity  /  since  the  impulsive  loading  may 
cause  the  cylinder  wall  to  experience  back-surface  spall.  Thus  in  this  case,  to  obtain 
the  correct  value  of  h  for  use  in  Equation  3.  we  would  first  need  to  evaluate  the 
functional  (or  algorithmic)  relationship  between  h  and  I .  BETAFACT  provides  for 
this  situation  by  allowing  the  user  to  write  another  FORTRAN  subroutine  (called 
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AUXALGO)  which  is  called  by  the  program.  This  subroutine,  too,  must  be  compiled 
and  linked  to  BETAFACT,  in  the  manner  described  in  Section  5,  to  become  available 
to  the  code.  Specific  details  concerning  the  structure  and  use  of  both  the  TF  and 
AUXALGO  subroutines  are  provided  in  Section  5. 

Once  the  executable  file  of  BETAFACT  containing  the  user’s  algorithm  is  pre¬ 
pared,  analysis  of  the  algorithm  by  the  code  proceeds  in  an  interactive  mode.  The 
program  prompts  the  user  for  all  additional  information  required  to  initially  analyze 
the  algorithm  or  re-analyze  it  according  to  user  specified  modifications.  The  addi¬ 
tional  information  required  include  identification  of  the  probability  distribution  type 
to  be  used  in  modeling  parameter  and  algorithm  uncertainties,  the  nominal  value 
for  each  parameter  with  associated  uncertainty  and  the  specification  (in  terms  of  K- 
factors)  of  that  uncertainty,  the  values  of  independent  parameters,  and  the  number 
of  times  the  algorithm  is  to  be  evaluated  to  develop  the  results  distribution. 

After  the  problem  to  be  processed  is  completely  defined,  BETAFACT  first  uses 
the  input  data  to  determine  the  precise  probability  distribution  for  each  parameter 
and  algorithm  with  uncertainty  which  will  be  used  to  model  the  uncertainty.  The 
code  then  proceeds  to  evaluate  the  algorithm  the  requested  number  of  times.  For 
each  evaluation,  a  different  normally  or  uniformly  (as  appropriate)  distributed  ran¬ 
dom  number  is  used  in  conjunction  with  a  parameter  nominal  value  and  associated 
probability  distribution  to  obtain  a  randomized  value  for  each  parameter  with  uncer¬ 
tainty.  Details  of  this  process  are  described  in  the  subsections  below.  The  algorithm 
is  then  evaluated  with  the  set  of  randomized  parameters,  uncertainty  is  applied  to 
the  result  of  the  computation  (if  required),  and  the  final  result  is  then  saved.  In  this 
manner  the  algorithm  results  distribution  is  generated. 

The  algorithm  results  distribution  contains  all  the  information  required  for  deter¬ 
mining  the  combined  uncertainty  associated  with  the  output  of  the  algorithm.  First 
the  results  distribution  is  sorted  and  ordered  from  lowest  to  highest  value.  Next  its 
mean,  variance,  and  standard  deviation  are  evaluated  using  Equations  13  and  14. 
These  quantities  (and  in  some  cases  the  low  and  high  extreme  values  of  the  algo¬ 
rithm  results  distribution)  are  all  the  information  needed  to  quantify  the  K-  factor 
uncertainty  associated  with  the  results  distribution.  These  K-factors  are  reported 
by  the  program  and  then  the  user  is  prompted  to  specify  whether  or  not  data  files 
containing  the  histogram  of  the  results  distribution  and  a  probability  distribution 
fit  to  the  histogram  are  to  be  generated.  These  files  contain  data  pairs  defining  the 
histogram  and  the  fit  and  can  be  plotted  with  software  external  to  BETAFACT.  The 
code  finally  prompts  the  user  to  interactively  modify  the  problem  and  re-analyze  it 
or  to  terminate  the  program  execution. 

The  above  description  of  the  algorithm  analysis  process  used  in  BETAFACT  em¬ 
phasizes  that  input  uncertainties  and  the  overall  algorithm  output  uncertainty  are 
in  terms  of  what  are  called  K-factors.  We  first  briefly  mentioned  the  K-factor  ap¬ 
proach  for  quantifying  uncertainty  in  our  description  of  the  lognormal  distribution  in 
subsection  2.2.2.  In  the  following  subsections,  we  will  see  that  the  meaning  of  the 
K-factor  specification  for  an  uncertainty  is  dependent  on  the  probability  distribution 
used  to  model  that  uncertainty.  It  is  crucial  that  the  user  of  BETAFACT  thoroughly 
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understand  the  K-factor  uncertainty  concept  in  order  to  be  able  to  properly  apply 
the  code  and  interpret  its  results. 

3.2  NORMALLY  DISTRIBUTED  UNCERTAINTY. 

Suppose  that  a  parameter  x  with  associated  uncertainty  is  known  to  be  or  is  modelled 
as  being  normally  distributed.  If  we  assume  the  mean  /x  and  standard  deviation  a 
of  this  parameter  distribution  are  known,  then  the  probability  distribution  for  the 
parameter  is  given  by  Equation  15.  We  recall  that  for  a  randomly  selected  value  of 
x,  95%  of  the  time  the  selected  value  will  fall  in  the  range  fx  —  1.96a  <  x  <  /x  +  1.96<r 
and  equally  as  often  above  /j  as  below  /i.  The  definition  of  the  K-factor  specification 
of  the  uncertainty  associated  with  this  parameter  is 

K  =  1.0  +  l.96—  (26) 

/* 

To  use  parameter  x  in  BETAFACT,  we  specify  it  as  being  normally  distributed 
and  enter  its  nominal  value  ANOM  =  fi  and  K-factor  uncertainty  AK  =  K.  (Here 
and  in  the  following  we  will  use  when  possible  variable  names  and  expressions  similar 
or  identical  to  those  actually  coded  in  BETAFACT.)  Then  the  low  and  high  side 
values  of  x  corresponding  to  the  95%  range  given  above  in  terms  of  /x  and  a  are, 
respectively,  ALO  —  2  *  ANOM  -  AH  I  and  AH  I  =  AK  *  ANOM.  Finally,  if  RN 
is  a  normally  distributed  random  number,  then  a  randomized  value  A  for  parameter 
x  is  given  by 

A  =  ANOM  +  RN*  {AH I  -  ANOM)/ 1.96  (27) 

Once  an  algorithm  results  distribution  is  generated  using  normally  distributed 
parameter  and/or  algorithm  uncertainties,  the  mean  /x  and  standard  deviation  a  of 
the  distribution  is  computed.  Equation  26  is  then  used  to  estimate  the  appropriate 
K-factor  uncertainty  which  should  be  associated  with  the  output  of  the  algorithm. 

3.3  LOGNORMALLY  DISTRIBUTED  UNCERTAINTY. 

Suppose  parameter  x  is  lognormally  distributed  with  mean  ANOM.  We  recall  from 
subsection  2.2.2  that  if  95%  of  the  time  a  randomly  selected  value  for  x  falls  in  the 
range  ANOM/ K  <  x  <  K  x  ANOM,  and  equally  as  often  below  ANOM  as  above 
it,  then  AK  =  K  is  the  K-factor  uncertainty  specification  for  x.  The  lognormal 
distribution  for  x  can  be  transformed  into  a  normal  distribution  in  y  =  lnx  with 
mean  and  standard  deviation  given  by 

/x  =  In  ANOM  (28) 

1.96 

The  low  and  high  side  values,  respectively,  of  x  corresponding  to  the  95%  range 
given  above  are  clearly  ALO  =  ANOM/ K  and  AHI  =  K  *  ANOM .  A  randomized 
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value  A  for  parameter  x  is  given  by 

A  =  ANOM*ALO**(RN/ 1.96)  (29) 

in  which  RN  is  a  normally  distributed  random  number. 

After  we  compute  the  mean  n  and  standard  deviation  a  of  an  algorithm  results 
distribution,  the  K-factor  estimate  to  be  associated  with  the  algorithm  output  is  given 
by  K  =  exp(1.96<r//i) 

3.4  BETA  DISTRIBUTED  UNCERTAINTY. 

Parameter  x  is  Beta  distributed  if  its  distribution  of  possible  values  is  given  by  Equa¬ 
tion  16.  This  distribution  is  not  as  conveniently  used  as  the  normal  and  lognor¬ 
mal  distributions.  In  fact,  we  don’t  use  the  Beta  probability  distribution  as  such  in 
BETAFACT,  but  instead  we  use  the  cumulative  probability  distribution  (CPD)  asso¬ 
ciated  with  a  given  Beta  distribution.  The  CPD  is  simply  the  total  area  (cumulative 
probability)  under  the  normalized  Beta  distribution  curve,  integrated  from  y  =  0  to 
y  =  y\  as  a  function  of  y' .  We  further  use  a  discrete  form  of  the  CPD  associated  with 
each  Beta  distribution  which  is  referenced  in  the  course  of  modeling  uncertainties. 
The  discrete  CPD  in  each  case  consists  of  a  list  of  the  21  values  for  y'  corresponding, 
respectively,  to  0%,  5%,  10%,  ...,  100%  of  the  cumulative  probability.  The  eleventh 
entry  in  such  a  list  corresponds  to  the  fraction  of  the  unit  interval  (0  to  1)  at  which 
50%  of  the  area  under  the  Beta  probability  distribution  curve  is  to  the  left  and  50% 
is  to  the  right  of  that  unit  interval  location. 

We  complicate  matters  further  by  offering  two  options  in  BETAFACT  regarding 
the  modelling  of  parameter  uncertainties  with  Beta  CPDs;  the  option  to  use  the  most 
appropriate  Beta  CPD  from  a  library  of  tabulated  distributions  or  to  have  the  code 
calculate  a  Beta  CPD  specifically  for  the  situation  at  hand.  The  former  option  is 
the  only  one  currently  offered  in  the  FAST  code.  The  second  option  increases  the 
accuracy  with  which  uncertainties  can  be  modelled  using  Beta  CPDs.  As  will  be  seen 
below,  however,  there  is  a  cost  associated  with  each  option. 

Before  considering  in  greater  detail  each  option,  we  note  that  both  require  the 
specification  of  two  K-factors  for  each  Beta  distributed  parameter.  A  K-factor  is 
required  to  characterize  each  side  or  the  distribution  relative  to  the  nominal  value 
ANOM  of  the  parameter.  The  low  side  K-factor  is  named  AKLO  in  BETAFACT 
and  the  high-side  K-factor  AKHI.  The  meaning  of  these  K-factors  is  that  100%  of 
the  possible  values  of  parameter  x  fall  in  the  range  ANOM/ AKLO  <  x  <  AK HI  * 
ANOM,  50%  of  the  time  below  ANOM  and  50%  of  the  time  above.  Thus  the 
low  and  high  side  extreme  values  of  the  parameter  are  ALO  =  ANOM / AK LO  and 
AH  I  =  AK  HI  x  ANOM,  respectively. 

A  randomized  value  A  for  the  parameter  is  obtained  in  a  somewhat  complicated 
process  using  the  21  discrete  y'  values  for  a  Beta  CPD  which  are  stored  in  an  array 
called  BD.  First  a  uniformly  distributed  random  number  RN  is  generated  and  used 
to  compute  the  number  F LI  =  20.0  *  RN  -I-  1.0.  This  latter  number  is  decomposed 
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into  its  integer  part  (0,1,2,. ..,20)  and  its  fractional  part.  The  integer  part  of  FLI  is 
used  to  select  2  adjacent  entries  of  the  BD  array.  For  instance,  if  the  integer  part  of 
FLI  is  z,  then  the  z’-th  and  (z  -F  1)  -th  BD  entries  are  selected.  The  fractional  part 
F  quantifies  a  fraction  of  the  distance  between  BD(z)  and  BD(z  +  1).  It  is  used  to 
calculate  the  value  of  variable  BDV; 

BDV  =  (l- F)*BD(i)  + F*BD(i+l)  (30) 

which  is  then  used  to  evaluate  the  randomized  value  for  the  variable  corresponding 
to  the  random  number  RN; 

A  =  ALO  +  BDV  *  (AH I  —  ALO)  .  (31) 


The  easiest  (and  less  accurate)  method  available  in  BETAFACT  for  modelling 
uncertainties  with  Beta  probability  distributions  is  to  use  the  Beta  CPDs  which  are 
tabulated  in  the  code.  There  are  81  of  these  CPDs.  Table  1  identifies  the  medians 
of  the  distributions  and  the  values  of  the  parameters  a  and  /?  which  were  used  to 
generate  the  distribution.  The  method  we  use  in  the  program  to  associate  a  specific 
tabulated  Beta  CPD  to  a  parameter  with  uncertainty  is  a  two  step  process.  First 
we  compute  the  median  (TEMP)  on  the  unit  interval  defined  by  the  nominal  value 
ANOM  and  K-factors  ( AKLO  and  AKHI)  of  the  uncertain  parameter: 


TEMP  = 


ANOM  -  ANOM/ AKLO  1  -  1  /AKLO 

AKHI  *  ANOM  -  ANOM  [AKLO  ~  AKHI  -  if  AKLO 


(32) 


We  then  determine  which  of  the  81  tabulated  Beta  CPDs  has  a  median  (11-th  entry) 
closest  to  TEMP.  The  closest  CPD  thereby  identified  is  used  to  model  the  uncertainty 
associated  with  the  parameter. 


The  more  accurate  method  in  BETAFACT  for  modeling  Beta  distributed  parame¬ 
ter  uncertainty  is  to  have  the  code  compute  the  Beta  CPD  which  corresponds  exactly 
to  the  distribution  desired.  This  is  achieved  by  specifying  for  an  uncertain  parame¬ 
ter  not  only  its  nominal  value  ANOM  and  associated  K-factors  AKLO  and  AKHI , 
but  also  the  standard  deviation  crq.  The  latter  must  satisfy  the  constraint  given  in 
Equation  21  which  BETAFACT  computes  and  reports  to  the  user.  With  these  four 
inputs,  then,  and  using  ALO  =  ANOM /AKLO,  AH  I  =  AKHI  x  ANOM ,  values 
for  m  and  a  are  computed  (Simons,  1988), 

ANOM  -ALO  <3q 

m  AH  I  -  ALO  a  AH  I  -  ALO  (  ^ 

and  then  Equation  20  is  evaluated  to  determine  corresponding  values  for  a  and  f3. 
Finally  Equation  16  (with  Equation  17)  is  evaluated  and  the  corresponding  Beta  CPD 
is  computed  and  stored  in  the  BD  array  used  for  modeling  the  parameter  uncertainty. 


BETAFACT  generates  three  estimates  for  the  K-factors  to  be  associated  with 
a  Beta  distributed  algorithm  results  distribution  with  mean  n  and  low  and  high 
extremes  X LO  and  XHI,  respectively.  These  are  the  low-side  K-factor  estimate, 
KLO  —  fi/XLO,  the  high-side  K-factor,  KHI  =  XHI/n,  and  the  average  of  the 
two  KAVE  =  0.5 {KLO  +  KHI). 
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3.5  UNIFORMLY  DISTRIBUTED  UNCERTAINTY. 

The  generalized  uniform  probability  distribution  for  parameter  z  is  given  by  Equa¬ 
tion  22.  This  distribution  is  characterized  in  terms  of  the  nominal  value  ANOM 
for  the  parameter  and  low  and  high  side  K-factors,  AKLO  and  AKHI,  respec¬ 
tively.  These  K-factors  mean  that  100%  of  the  possible  range  for  parameter  z  is  given 
by  (ANOM/ AKLO  <  x  <  AKHI  *  ANOM)  with  50%  of  the  distribution  below 
ANOM  and  50%  above.  As  in  the  case  of  a  Beta  distributed  parameter,  the  extreme 
values  for  the  parameter  are  ALO  =  ANOM  /  AKLO  and  AH  I  =  AKHI  *  ANOM. 

A  uniformly  distributed  random  number  RN  is  used  to  obtain  a  randomized  value 
A  for  parameter  z.  If  the  value  of  RN  is  less  them  or  equal  to  0.5,  then 

A  =  ALO  +  2.0  *  RN  *  ( ANOM  -  ALO)  (34) 

else 

A  =  ANOM  +  2.0  *  (RN  -  0.5)  *  (AH I  -  ANOM)  (35) 

The  same  three  K-factor  estimates  produced  by  BETAFACT  for  Beta  distributed 
algorithm  results  are  also  generated  for  algorithm  results  distributions  computed  using 
uniform  distributions  to  model  parameter  and  algorithm  uncertainties. 
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Table  1.  Tabulated  Beta  cumulative  probability  distribution  means  and  a  and  f3 
parameters. 
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SECTION  4 

REPRESENTATIVE  ANALYSES 


Besides  having  been  used  extensively  in  our  lethality  assessment  activities  for  the 
LTH-3  Program,  BETAFACT  has  been  exercised  on  many  simple  problems  to  verify 
the  implementation  in  the  code  of  the  probability  distributions  described  in  Section  2 
and  the  methodology  of  Section  3.  It  has  also  been  used  to  analyze  more  complicated 
problems  which  more  completely  illustrate  the  several  uses  of  the  code.  The  results  of 
several  simple  BETAFACT  analyses  are  reviewed  in  subsection  4.1.  The  application  of 
the  code  in  the  analysis  of  a  more  realistic  algorithm  is  the  central  theme  of  subsection 
4.2. 


4.1  VERIFICATION  ANALYSES. 

We  consider  the  simplest  possible  algorithm  r  =  x  in  analyses  to  verify  proper  im¬ 
plementation  of  the  theory  and  methodology  in  BETAFACT.  This  algorithm  has  the 
output  r  and  the  single  input  parameter  x.  The  input  parameter  is  assumed  to  have 
associated  uncertainty.  Towards  the  end  of  this  subsection,  we  also  consider  the  case 
in  which  the  algorithm  has  associated  uncertainty. 

A  listing  of  subroutine  TF  used  to  input  this  simple  algorithm  to  BETAFACT  is 
given  in  Appendix  A.  The  subroutine  locally  references  variable  x  (which  is  passed  to 
the  subroutine  as  the  first  element  A(l)  in  array  A  as  ADUM.  The  subroutine  also 
uses  the  independent  parameter  XV  AL(1)  which  has  an  interactively  specified  value 
of  1.0  and  hence  has  no  effect  on  the  algorithm  output  r  =  YV AL.  More  details 
about  formulating  subroutine  TF  axe  given  in  subsection  4.2. 

In  the  situation  when  uncertainty  is  associated  only  with  the  input  parameter  x, 
we  know  that  the  distribution  of  the  algorithm  results  should  be  identical  to  that  used 
to  model  the  input  parameter  uncertainty.  In  the  following  the  latter  is  modelled,  in 
turn,  using  each  of  the  four  probability  distribution  types  available  in  BETAFACT. 
We  assume  a  nominal  value  for  the  parameter  of  50,000  and  analyze  the  problem  for 
a  few  different  values  of  the  K-factor  uncertainty.  The  algorithm  is  evaluated  5000 
times  in  each  of  the  referenced  analyses  in  order  to  generate  the  results  distribution. 
A  fewer  number  of  algorithm  evaluations  (e.g.,  1000)  gives  results  not  substantially 
different  from  those  we  present  below. 

The  initial  set  of  analysis  results  we  consider  are  obtained  when  the  parameter 
uncertainty  is  modelled  as  normally  distributed  with  a  K-factor  specification  of  1.25 
or  1.50.  Plots  of  the  computed  (solid  line)  and  fitted  (dash  line)  results  distribution 
for  the  two  K-factor  specifications  are  shown  in  Figure  5.  The  figures  also  report  the 
random  number  generator  seed  integer  used  at  the  start  of  each  analysis.  The  value 
of  this  seed  integer  is  purposely  and  arbitrarily  varied  from  one  analysis  to  another. 

The  plots  in  Figure  5  show  that  the  computed  results  distribution  and  the  fitted 
distribution  generally  have  similar  shapes.  No  significant  effort  is  made  in  BETAFACT 
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Figure  5.  Simple  algorithm  results  distributions  for  K  =  1.25  (top)  and  K  =  1.50 
(bottom)  and  normally  distributed  uncertainty. 
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to  obtain  the  “best  fit”  to  the  computed  results  distribution;  a  simple  scaling  proce¬ 
dure  is  employed  to  approximately  match  the  peak  of  the  distribution.  The  purpose 
of  the  plotted  distributions  is  more  of  a  qualitative  rather  than  a  quantitative  check 
of  the  results.  What  we  are  really  after  in  using  BETAFACT  is  the  K-factor  as¬ 
sociated  with  the  results  distribution  and  the  mean  of  the  distribution.  These  two, 
together  with  the  probability  distribution  type,  are  sufficient  to  completely  character¬ 
ize  the  uncertainty  distribution  associated  with  the  algorithm  output  for  subsequent 
use  in  FAST.  Thus  to  verify  that  the  code  is  functioning  properly,  we  must  compare 
the  computed  statistics  of  the  results  distribution  to  the  corresponding  theoretical 
values. 

In  defining  this  simple  problem,  we  specified  the  nominal  value  of  the  parameter 
as  y.  =  50,000  and  its  K-factor  uncertainty  as  1.25  or  1.50.  Equation  26  enables 
us  to  relate  these  values  to  the  corresponding  standard  deviation  a  of  the  normal 
probability  distribution  which  exactly  models  the  parameter  and  hence  the  algorithm 
results  uncertainty.  Doing  so,  we  find  that  the  exact  values  for  a  are  6378  for  K  =  1.25 
and  12,755  for  K  =  1.50.  The  values  actually  computed  by  BETAFACT  for  p  and 
a  and  for  the  K-factor  estimate  of  the  results  distribution  are  reported  in  Table  2. 
Compared  to  the  theoretical  values,  the  computed  results  are  seen  to  be  quite  good. 
Two  reasons  the  comparisons  aren’t  better  are  that  an  infinite  number  of  algorithm 
evaluations  were  not  used  to  generate  the  results  distribution  (a  small  effect  here)  and 
only  100  distinct  frequency  of  occurrence  bins  were  used  to  construct  the  histograms 
(most  significant  effect). 

The  above  simple  analysis  problem  was  repeated,  this  time  modelling  the  param¬ 
eter  uncertainty  distribution  as  lognormal.  The  plotted  results  from  these  analyses 
are  shown  in  Figure  6.  The  equation  at  the  end  of  subsection  3.3  enables  us  to  relate 
the  specified  K-factor  and  mean  to  the  corresponding  theoretical  standard  deviation. 
The  theoretical  values  are  compared  to  the  BETAFACT  computed  values  in  Table 
2.  The  latter  are  within  a  few  percent  of  the  theoretical  values.  The  reasons  given 
above  for  the  normal  distribution  computed  results  discrepancy  apply  here  as  well. 

We  next  analyzed  the  simple  algorithm  using  Beta  cumulative  probability  distri¬ 
butions  to  model  the  parameter  uncertainty.  The  two  analyses  described  above  were 
each  performed  twice,  first  using  tabulated  Beta  CPDs  and  then  using  BETAFACT 
computed  distributions.  The  reader  may  recall  that  a  standard  deviation  for  the 
uncertain  parameter  must  bo  specified  to  allow  the  use  of  computed  rather  than  tab¬ 
ulated  distributions.  The  standard  deviations  computed  when  the  uncertainties  were 
modelled  with  tabulated  distributions  were  used  as  the  exact  input  standard  devia¬ 
tions  for  the  calculated  distribution  analyses.  The  computed  results  for  both  sets  of 
analyses  are  presented  in  Table  2.  Figure  7  shows  the  distribution  plots  obtained  in 
the  tvo  computed  Beta  distribution  analyses.  The  tabulated  distribution  results  are 
nearly  identical  to  those  shown  in  the  figure. 

Finally,  we  analyzed  the  simple  algorithm  example  using  a  uniformly  distributed 
uncertainty  model.  Three  combinations  of  K-factors  were  considered.  For  the  first 
two  K-factor  combinations  ((1.25,1.20)  and  (2.00,1.50))  Equation  25  provides  the 
theoretical  value  of  the  corresponding  standard  deviation.  Equation  24  must  be  used 
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Table  2.  Simple  algorithm  analysis  results. 


Distribution 

Specified 

K- Factors 

Exact 

Computed 

Type 

KLO 

KHI 

<T 

a 

KLO 

KHI 

Normal 

1.25 

- 

6378 

49976 

6249 

1.245 

- 

1.50 

- 

12755 

50152 

12769 

1.499 

- 

Lognormal 

1.25 

- 

5692 

50316 

5829 

1.255 

- 

1.50 

- 

10344 

51167 

10674 

1.505 

Beta* 

-Tabulated 

1.25 

1.25 

50169 

4223 

1.254 

1.245 

1.50 

1.50 

- 

50434 

7877 

1.513 

1.487 

-Calculated 

1.25 

1.25 

4223 

50054 

4489 

1.251 

1.249 

1.50 

1.50 

7877 

50111 

8284 

1.503 

1.496 

Uniform 

1.25 

1.20 

5774 

49998 

5806 

1.250 

1.200 

2.00 

1.50 

14434 

50002 

14500 

1.999 

1.500 

2.00 

2.00 

21949 

48995 

23046 

1.960 

2  .041 

f-  Analyses  performed  first  with  tabulated  (unspecified  cr)  and  then  with 
calculated  (specified  a)  Beta  cumulative  probability  distributions. 
Theoretical  mean  is  50000  for  all  cases. 
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Simple  algorithm  results  distributions  for  K  =  1.25  (top)  and  A 
(bottom)  and  Beta  distributed  uncertainty. 


to  evaluate  the  theoretical  standard  deviation  associated  with  the  third  K-factor  pair 
(2.00, 2.00).  The  theoretical  and  BETAFACT  computed  values  for  the  mean,  standard 
deviation,  and  K-factors  of  these  analyses  axe  given  in  Table  2.  The  agreement 
between  computed  and  theoretical  values  is  quite  good.  Figures  8  and  9  show  plots 
of  the  computed  and  fitted  results  distributions. 

While  the  plots  in  Figures  5  through  9  are  useful  in  visualizing  the  output  from 
BETAFACT,  it  is  the  comparisons  evident  in  Table  2  which  indicate  that  the  code  is 
functioning  properly  in  the  analyses  of  the  simple  algorithm  r  =  x.  These  comparisons 
show  that  the  code  adequately  replicates  the  K-factors  and  nominal  value  specified 
for  the  uncertain  algorithm  parameter.  It  is  somewhat  less  successful  in  reproducing 
the  theoretical  standard  deviation.  The  limited  structure  (100  bins)  of  the  histograms 
compiled  from  the  computed  algorithm  results  distribution  is  believed  to  be  the  main 
source  for  the  error  in  the  computed  standard  deviations. 

As  a  final  verification  analysis,  we  considered  the  case  in  which  both  the  input 
parameter  and  the  algorithm  have  associated  uncertainty.  If  both  of  these  uncertainty 
sources  are  modeled  as  normally  distributed,  then  the  theoretical  result  is  that  the 
variance  (standard  deviation  squared)  of  the  results  distribution  is  equal  to  the  sum 
of  the  variance  associated  with  the  input  parameter  and  that  associated  with  the 
algorithm.  Specifying  the  parameter  nominal  value  again  as  50, 000,  the  parameter 
K-factor  as  1.25,  and  the  algorithm  K-factor  also  as  1.25,  the  parameter  standard 
deviation  is  6378  (as  is  that  of  the  algorithm),  and  the  results  distribution  has  a 
theoretical  standard  deviation  of  9020  (and  corresponding  K-factor  of  1.354).  When 
this  problem  is  analyzed  with  BETAFACT  using  5000  evaluations,  we  obtain  the 
following  computed  values;  p.  —  50, 067,  a  =  9088,  and  K  =  1.356.  The  computed 
values  are  in  very  good  agreement  with  the  theoretical  results.  Plots  of  the  computed 
histograms  for  this  case  are  shown  in  Figure  10. 

4.2  REALISTIC  ALGORITHM  EXAMPLE. 

The  analysis  results  presented  in  the  previous  subsection  give  us  confidence  that  the 
probabilistic  theory  we  are  using  is  correctly  implemented  in  BETAFACT.  In  the 
current  subsection  we  illustrate  how  the  code  is  used  to  analyze  a  realistic  algorithm. 
Specifically  we  consider  a  slightly  modified  version  of  the  algorithm  presented  previ¬ 
ously  in  Equation  3.  The  modified  version  of  this  basic  algorithm  is  as  follows: 

£._S_ 

2perhl'5 

This  algorithm  provides  an  estimate  for  the  peak  strain  £  in  a  metallic  cylinder  in 
terms  of  the  parameters  on  the  right  hand  side  of  the  equation.  These  parameters 
are  material  density  p,  material  ultimate  strength  <r,  thickness  hr  and  peak  impulse 
intensity  Ir.  We  assume  that  each  parameter  used  in  the  right  hand  side  of  the 
algorithm,  as  well  as  the  algorithm  itself,  has  associated  uncertainty. 

The  algorithm  given  in  Equation  36  is  different  from  the  Equation  3  algorithm 
in  three  respects.  First,  we  have  arbitrarily  changed  the  exponent  of  the  thickness 


(36) 
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Number  of  Occurrences  Number  of  Occurrences 


Number  of  Occurrences 


parameter  from  2  to  1.5.  We  have  then  added  the  subscript  r  to  h  and  7  to  indicate 
that  the  actual  thickness  of  the  cylinder  wall  and  the  actual  impulse  intensity  which 
loads  the  structure  are  not  used  in  the  algorithm.  Instead  we  assume  we  must  first 
compute  these  parameter  values  from  auxiliary  algorithms  before  the  proper  values 
can  be  used  in  the  strain  algorithm.  In  this  example  analysis,  we  assume  the  following 
forms  for  these  auxiliary  algorithms; 

Ir  =  [3.5  +  0.08(7  -  3.5)]  (2h)1/2  (37) 

and 

hT  =  h-  [0.05  +  0.001(7  -  3.5)]  (2  h)1/2  (38) 

In  these  equations  h  is  the  actual  thickness  of  the  cylinder  wall  before  the  impulsive 
load  is  applied  and  7  is  the  peak  intensity  of  the  actual  applied  impulsive  load.  The 
two  auxiliary  algorithms  are  used  here  to  model  the  situation  in  which  the  applied 
impulsive  load  causes  the  inner  surface  of  the  cylinder  wall  to  spall  if  the  impulse 
intensity  is  sufficiently  high.  When  spall  occurs,  the  thickness  of  the  residual  cylinder 
wall  is  reduced  and  some  of  the  momentum  originally  delivered  to  the  wall  is  car¬ 
ried  away  by  the  spall.  The  auxiliary  algorithms  enable  us  to  estimate  the  residual 
thickness  of  the  wall,  if  spall  occurs,  and  the  momentum  remaining  in  it.  These  are 
the  required  inputs  for  the  algorithm  of  Equation  36.  The  auxiliary  algorithms  are 
intended  to  be  used  only  if  the  actual  applied  impulse  intensity  is  greater  than  3.5. 

Summarizing  the  problem  defined  so  far,  we  have  a  basic  algorithm  which  has  four 
input  parameters  with  associated  uncertainty  .  The  algorithm  itself  is  also  assumed 
to  have  associated  uncertainty.  Two  of  the  parameters  used  in  the  algorithm  must 
first  be  computed  from  auxiliary  algorithms.  We  now  describe  in  some  detail  how  all 
these  algorithms  are  made  available  to  the  code. 

The  basic  algorithm,  Equation  36,  becomes  available  to  the  code  via  the  user  sup¬ 
plied  subroutine  TF.  In  order  to  prepare  the  required  subroutine,  we  first  must  decide 
in  what  order  we  wish  to  refer  to  the  uncertain  parameters.  This  is  necessary  since 
the  randomized  value  for  each  parameter  used  in  each  evaluation  of  the  algorithm  is 
stored  in  a  particular  location  in  an  array  called  A.  To  code  the  basic  algorithm  in  the 
present  example  analysis,  we  elect  to  associate  the  randomized  value  for  p  with  array 
element  A(l),  the  randomized  value  for  a  with  A{ 2),  7r  with  A(3),  and  hT  with  >4(4). 
To  accomplish  this  association,  we  simply  input  the  nominal  values  (and  associated 
K-factors)  for  these  parameters  in  this  order  when  prompted  for  them  during  pro¬ 
gram  execution.  Since  both  7r  and  hT  will  be  computed  by  the  auxiliary  algorithms, 
the  nominal  values  we  input  for  these  two  parameters  can  be  anything  (provided  we 
don’t  intend  to  use  the  input  nominals  in  the  auxiliary  algorithm  subroutine).  We 
typically  input  unity  (1.0)  as  the  nominal  value  for  parameters  such  as  7r  and  hT. 
The  actual  nominal  values  for  p  and  a  are  input.  All  of  these  nominal  values  are 
stored  by  BETAFACT,  in  the  order  of  their  entry,  in  array  ANOM.  During  program 
execution,  ANOM(l)  is  used  in  the  computation  of  A(l)  (the  randomized  parameter 
value),  ANOM{2)  for  A(2),  etc. 

If  the  basic  algorithm  under  consideration  is  assumed  not  to  have  associated  un¬ 
certainty,  then  the  nominal  values  described  above  are  all  the  nominals  expected  by 
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the  code.  With  specified  algorithmic  uncertainty,  however,  it  is  necessary  to  enter 
an  additional  nominal  value  (with  a  value  which  must  be  1.0  in  this  case)  and  the 
associated  K-factor  specifications.  In  our  example,  this  nominal  value  can  be  thought 
of  as  another  parameter  which  is  assigned  to  array  element  A(5)  (or  ANOM( 5)).  If 
algorithmic  uncertainty  is  to  be  applied  in  a  BETAFACT  analysis,  its  nominal  value 
(1.0)  and  associated  K-factors  must  be  the  last  set  of  uncertain  parameter  nominal 
inputs  to  the  code. 

In  most  instances  (though  not  in  the  present  case)  we  use  the  values  of  independent 
parameters  in  the  evaluation  of  the  basic  algorithm  in  subroutine  TF.  The  values 
of  these  independent  parameters  are  contained  in  array  XV AL.  The  order  of  the 
independent  parameter  values  in  this  array  is  determined  by  the  order  in  which  the 
independent  parameter  values  axe  input  at  the  appropriate  prompt  by  the  code.  In 
the  present  example  there  are  two  independent  parameters;  the  actual  applied  impulse 
intensity  I  which  we  choose  to  enter  first  and  which  is  thus  assigned  to  XV AL(  1)  and 
the  pre-spall  thickness  h  which  is  stored  in  XV AL{ 2).  The  entries  in  array  XV AL 
axe  not  changed  during  the  execution  of  BETAFACT  unless  the  user  selects  an  option 
which  enables  them  to  be  altered  interactively. 

Once  we  have  picked  an  ordering  of  the  uncertain  parameters  and  the  independent 
parameters,  we  can  then  proceed  to  formulate  subroutine  TF.  The  listing  of  TF  which 
codes  the  basic  algorithm  of  this  example  problem  is  provided  in  Appendix  A.  The 
algorithm  is  simple  and  its  coding  is  straightforward.  To  facilitate  checking  of  the 
algorithm  coding,  we  use  local  variable  names  in  the  subroutine  mnemonic  of  the 
respective  variables  in  the  algorithm.  For  instance,  we  set  ARHO  =  A(l)  since  array 
element  A(l)  is  the  randomized  value  of  the  density,  ASIG  =  A( 2)  which  is  the 
randomized  ultimate  stress,  etc.  (In  the  listing  given  for  TF  we  also  use  a  constant 
multiplier  1.0  x  106.  This  multiplier  serves  to  convert  impulse  intensity  /,  input  in 
ktaps,  to  taps  for  consistency  with  the  units  of  the  other  algorithm  input  parameters.) 

Subroutine  AUXALGO,  which  contains  the  auxiliary  algorithms,  is  prepared  in 
much  the  same  way  as  subroutine  TF.  A  listing  of  the  AUXALGO  subroutine  used  in 
the  present  analysis  is  given  in  Appendix  A.  It  is  important  to  note  that  we  use  the 
nominal  parameter  value  array  ANOM  and  not  the  array  A  (randomized  parameter 
value  array)  in  order  to  compute  the  auxiliary  algorithms.  Once  again  mnemonic  local 
variable  names  are  recommended  to  help  the  user  understand  and  verify  the  coding 
of  the  auxiliary  algorithms.  For  instance,  in  this  subroutine  we  use  AI  =  XV AL(  1) 
to  locally  represent  the  applied  impulse  intensity  and  AH  =  XV AL( 2)  the  original 
cylinder  wall  thickness. 

After  both  subroutine  TF  and  AUXALGO  are  prepared,  we  compile  and  link 
them  to  the  BETAFACT  object  file  to  generate  the  executable  version  of  the  code 
which  applies  specifically  to  our  problem.  The  actual  VAX/ VMS  commands  needed 
for  compiling  and  linking  the  subroutines  are  described  in  Section  5.  We  can  then 
proceed  to  run  the  code  interactively. 

The  first  example  analysis  we  consider  of  the  realistic  algorithm  described  above 
is  determination  of  the  overall  uncertainty  associated  with  the  output  of  the  algo- 
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Figure  11.  Realistic  algorithm  results  distribution  for  original  K-factor  specifications 
and  lognormally  distributed  uncertainty. 

rithm.  As  specific  input  we  use  p  =  2.7  gm/cm3  with  Kp  =  1.05,  a  =  3.45  x  109 
dyne/cm2(50ksi)  with  Ka  —  1.10,  IT  (nominal)  =  1.0  with  Ki  —  1.25,  hr  (nominal) 
=  1.0  with  Kh  =  1.40,  and  algorithm  (nominal)  =  1.0  with  Kaig0  =  1-25.  We  choose 
to  model  all  the  uncertainties  with  lognormal  distributions.  For  the  independent 
parameters  we  specify  I  (applied)  =  10.0  ktap  and  ^(original)  =  0.38cm.  We  note 
that  the  specifications  for  K[  and  Kh  are  the  K-factors  associated  with  the  outputs 
of  the  auxiliary  algorithms,  IT  and  hr,  respectively.  These  K-factors  represent  the 
effective  uncertainty  associated  with  the  auxiliary  algorithms,  including  all  impor¬ 
tant  contributions  to  the  overall  algorithmic  uncertainty.  Two  previous  BETAFACT 
analyses,  one  for  each  auxiliary  algorithm,  axe  assumed  to  have  been  used  to  generate 
the  combined  uncertainty  specification  for  each  auxiliary  algorithm. 

Using  the  values  described  above,  evaluating  the  algorithm  5000  times  to  generate 
the  results  distribution,  and  then  computing  the  statistics  of  the  distribution,  we  find 
that  its  mean  is  p  =  3.69  x  10~ 3 ,  its  standard  deviation  is  a  =  1.42  x  10~3  and 
the  associated  K-factor  uncertainty  is  Ke  =  2.13.  A  plot  of  the  results  distribution 
(continuous  line  curve)  and  fit  to  it  (dash  line  curve)  is  shown  in  Figure  11.  The 
theoretical  mean  for  this  problem  is  3.47  x  10~ 3,  which  is  about  6%  less  than  the 
computed  mean. 

A  very  useful  feature  of  BETAFACT  is  that,  once  a  problem  has  been  fully  defined 
for  a  set  of  input  parameters  and  associated  K-factors,  it  is  easy  to  re-analyze  the 
problem  using  varied  parameter  values  and/or  uncertainties.  Thus  we  can  use  the 
code  to  perform  algorithm  sensitivity  studies. 
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As  a  first  example  of  such  studies,  we  determine  which  source  of  uncertainty  in 
our  problem  is  the  most  significant  contributor  to  the  overall  combined  uncertainty. 
Clearly  the  most  significant  source  is  either  the  uncertainty  associated  with  the  basic 
algorithm  or  with  one  of  the  two  auxiliary  algorithms.  If  we  re-analyze  the  prob¬ 
lem  described  above,  but  with  Ki  =  1.01  (ie.,  very  little  uncertainty),  we  find  that 
the  code  estimates  the  combined  uncertainty  to  be  K  —  1.77;  17%  less  than  Ke 
given  above.  For  Kh  =  1.01  and  all  other  inputs  as  originally  specified,  we  find  that 
K  =  1.68;  21%  less  than  Ke.  Finally,  for  Kaigo  =  1.01  and  all  other  inputs  with  their 
original  values,  we  compute  K  =  2.02;  only  5%  less  than  Ke.  Thus  the  uncertainty 
associated  with  the  auxiliary  algorithm  which  computes  the  cylinder  wall  thickness 
is  the  most  significant,  that  associated  with  the  residual  impulse  intensity  algorithm 
only  slightly  less  significant,  and  that  attributed  to  the  basic  algorithm  in  relative 
terms  almost  insignificant.  While  this  ordering  of  uncertainty  significance  may  have 
been  anticipated  since  Kh  is  numerically  the  largest  of  the  three  uncertainties  con¬ 
sidered,  the  greater  exponent  (2.0)  on  Ir  in  the  basic  algorithm  compared  to  that 
of  hr  (1.5)  nearly  reverses  the  ordering.  A  more  complicated  algorithm  would  make 
determination  by  inspection  of  the  most  significant  uncertainty  source  a  much  more 
difficult  task.  However,  it  would  be  straightforward  with  BETAFACT. 

Now  suppose  we  are  asked  which  of  the  following  is  more  beneficial  in  reducing 
the  overall  combined  uncertainty;  reducing  Kj  from  1.25  (its  original  value)  to  1.10 
or  reducing  Kh  from  1.40  (its  original  value)  to  1.25?  Further  suppose  the  reductions 
will  be  achieved  via  test  programs  with  the  former  estimated  to  cost  S200K  and  the 
latter  $100K.  Which  test  program  should  be  pursued?  BETAFACT  enables  us  to 
quantify  justifiable  answers  to  these  questions.  The  code  can  thus  be  employed  as  a 
useful  management  and  decision  making  tool. 

First  consider  the  overall  uncertainty  reduction  obtained  by  reducing  the  con¬ 
tributing  uncertainties.  Re-analyzing  the  original  problem  with  Kj  reduced  from 
1.25  to  1.10,  we  find  that  the  combined  K-factor  decreases  from  2.13  to  1.83,  a  14% 
reduction.  If  instead  we  reduce  Kh  from  1.40  to  1.25,  we  find  that  Ke  decreases  from 
2.13  to  1.88,  a  12%  reduction.  Thus  the  maximum  benefit  is  obtained  by  reducing 
K[  in  this  case. 

Now  consider  the  cost  factor.  On  the  one  hand,  reducing  Ke  by  reducing  Kj  costs, 
on  the  average,  $14. 3K  per  percentage  point.  On  the  other,  reducing  Ke  by  improv¬ 
ing  Kh  costs  only  S8.3K,  on  the  average,  per  percentage  point.  Thus  a  percentage 
point  improvement  in  the  latter  case  costs  only  58%  what  it  does,  on  the  average, 
in  the  former.  From  the  viewpoint  of  maximized  uncertainty  reduction,  we  would 
recommend  the  option  which  reduces  K[.  However,  from  a  cost  viewpoint,  clearly  it 
would  be  more  economical  to  pursue  the  Kh  reduction  option. 

Other  sensitivity  studies  which  can  be  accomplished  using  BETAFACT  are  prob¬ 
ably  evident  to  the  reader.  The  above  simple  example  clearly  demonstrates  that  the 
code  is  adept  at  performing  such  quantitative  analysis. 
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SECTION  5 

BETAFACT  USER’S  MANUAL 


The  analysis  of  an  algorithm  with  BETAFACT  is  a  two  step  process;  preparation 
of  the  problem  definition  (including  formulation  of  user  supplied  subroutines)  and 
analysis  of  the  problem  with  the  code.  Since  BETAFACT  executes  interactively,  the 
bulk  of  the  effort  involved  in  analyzing  a  new  algorithm  is  incurred  in  defining  the 
problem. 

In  the  following,  we  describe  some  preliminaries  to  the  execution  of  the  code  (sub¬ 
section  5.1),  review  the  steps  involved  in  formulating  a  problem  definition  (subsection 
5.2)  and  describe  in  detail  the  interactive  execution  of  BETAFACT  (subsection  5.3). 
Many  details  of  the  problem  definition  phase  are  described  more  completely  in  pre¬ 
vious  sections  of  this  report  (especially  Sections  3  and  4)  and  will  be  repeated  only 
briefly  here.  The  reader  should  refer  to  the  previous  report  sections  for  the  additional 
details. 

5.1  PRELIMINARIES  TO  EXECUTION. 

We  first  assume  the  user  has  installed  the  file  containing  the  BETAFACT  source  code 
(named  for  reference  here  as  BETAFACT. FOR)  and  the  object  file  (BETAFACT. OBJ) 
on  the  computer  to  be  used  in  the  algorithm  analyses.  Installation  here  means  simply 
having  a  copy  of  each  file  in  the  user’s  local  file  area.  A  listing  of  the  source  code  is 
provided  as  Appendix  B  to  this  report.  It  can  be  used  to  verify  the  completeness  of 
the  user’s  source  file. 

If  the  user  does  not  have  file  BETAFACT. OBJ  available,  but  does  have  the  source 
file,  an  object  file  can  be  generated  with  the  command 

FOR/CONTINUATIONS=99  BETAFACT.FOR 

The  above  command  applies  to  compilation  of  the  code  on  VAX  type  computers 
operating  under  VMS.  Since  BETAFACT  is  coded  using  standard  FORTRAN-77, 
probably  it  can  be  readily  compiled  on  any  computer  (such  as  a  PC)  with  a  FOR¬ 
TRAN  compiler.  A  command  similar  to  the  above  would  have  to  be  issued  on  the 
non- VAX  computer  to  generate  the  desired  BETAFACT  object  file.  The  CONTINU¬ 
ATIONS  parameter  in  the  above  command  is  necessary  since  some  data  statements 
in  the  BETAFACT  source  continue  over  50  or  more  lines  of  coding. 

In  order  to  obtain  an  executable  version  of  BETAFACT  which  is  specifically  for 
analysis  of  the  user’s  algorithm,  the  algorithm  must  be  made  available  to  the  code 
via  the  user  supplied  subroutine  named  TF  (which  stands  for  Transfer  Function  in 
the  FAST  terminology).  If  the  basic  algorithm  under  consideration  happens  to  use 
auxiliary  algorithms  (i.e.,  algorithms  which  compute  parameter  values  which  are  part 
of  the  input  to  the  basic  algorithm),  these  must  be  made  available  to  the  code  via  the 
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user  supplied  subroutine  named  AUXALGO  (for  auxiliary  algorithm).  These  two  sub¬ 
routines  contain  the  FORTRAN  coding  of  the  user’s  basic  and  auxiliary  algorithms, 
respectively. 

Once  the  user’s  subroutines  are  written  and  compiled,  they  must  be  linked  to  the 
BETAFACT  object  file  as  follows: 

1.  If  auxiliary  algorithms  are  used,  then  link  with 


LINK  BETAFACT, TF, AUXALGO 
2.  If  auxiliary  algorithms  are  not  used,  then  simply  link  with 


LINK  BETAFACT, TF 

A  nonfatal  error  message  will  be  obtained  with  the  second  LINK  above.  This 
merely  notifies  the  user  that  AUXALGO  is  referenced  in  BETAFACT  but  is  not 
defined  since  no  subroutine  AUXALGO  was  linked.  The  code  will  execute  properly 
even  though  this  nonfatal  error  was  encountered  during  linking. 

The  result  of  the  above  linking  procedures  is  a  file  named  BETAFACT.EXE  which 
is  the  executable  version  of  the  code  tailored  specifically  for  analysis  of  the  user’s 
algorithm.  To  begin  interactive  execution  of  this  program  file,  it  is  only  necessaxy  to 
issue  the  command  (on  VAX  type  computers). 

RUN  BETAFACT 

From  this  point  on,  the  user  need  only  respond  to  the  prompts  issued  by  the  code. 

5.2  PROBLEM  DEFINITION. 

The  user  needs  to  understand  three  single  index  arrays  and  one  variable  used  by 
BETAFACT  in  order  to  write  useful  user  subroutines  for  the  code.  The  three  arrays 
are  named  ANOM(i),  XV AL(i)  and  A(i)  and  the  variable  YV AL.  The  meaning  of 
each  in  the  code  is  as  follows: 

1.  ANOM(i)  -  This  array  contains  the  nominal  values  of  the  parameters  with 
associated  uncertainty  which  are  used  in  the  basic  algorithm.  Each  source  of 
uncertainty  is  identified  with  a  specific  entry  in  array  ANOM.  The  user  defines 
the  association  of  an  uncertain  parameter  with  say  the  z-th  entry  in  ANOM  by 
entering  it  as  the  z-th  parameter  value  when  the  code  prompts  for  nominal  values 
of  parameters  with  associated  uncertainty.  If  the  actual  parameter  nominal 
value  is  to  be  determined  by  an  auxiliary  algorithm  in  subroutine  AUXALGO 
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(see  subsection  4.2  for  a  specific  example  of  this  situation),  the  nominal  value 
specified  interactively  for  the  parameter  can  be  any  value,  provided  the  input 
value  isn’t  used  in  AUXALGO.  We  recommend  using  a  nominal  value  of  1.0 
for  such  parameters.  If  the  user’s  basic  algorithm  is  to  have  an  associated 
uncertainty  (again  see  subsection  4.2  for  an  example),  this  constitutes  a  source 
of  uncertainty  which  is  treated  in  BETAFACT  as  a  parameter  with  a  nominal 
value  which  must  be  1.0.  This  parameter  nominal  must  also  be  the  last  uncertain 
parameter  nominal  value  entered  interactively  by  the  user. 

2.  XV AL(i)  -  This  array  contains  the  values  of  parameters  used  in  the  user’s 
algorithms  which  do  not  have  associated  uncertainty.  They  are  also  referred 
to  as  independent  variables  or  parameters.  The  value  of  the  i-th  independent 
parameter  is  associated  with  the  i-th  entry  of  array  XV AL  by  it  being  the  2-th 
value  entered  at  the  BETAFACT  prompt  for  the  nominal  values  of  independent 
parameters. 

3.  A(i)  -  This  array  contains  in  its  i-th  location  the  current  randomized  value  of 
an  uncertain  parameter  corresponding  to  the  nominal  parameter  value  stored 
as  the  i-th  entry  in  array  ANOM.  Section  3  describes  how  these  randomized 
values  are  generated. 

4.  YV AL  -  This  variable  has  the  value  computed  for  the  output  of  the  user’s  basic 
algorithm  for  the  current  set  of  randomized  input  parameters  but  before  basic 
algorithm  uncertainty  has  been  applied. 

The  arrays  ANOM  and  XV AL  are  accessible  in  subroutine  AUXALGO.  Indeed, 
the  purpose  of  AUXALGO  is  to  assign  one  or  more  of  the  values  of  array  ANOM. 
All  three  arrays  can  be  used  in  subroutine  TF.  The  user  must  be  careful  not  to  change 
any  of  the  array  entries  (particularly  those  of  ANOM  and  XV AL)  in  that  subroutine, 
however.  The  value  of  the  evaluated  algorithm  returned  by  subroutine  TF  must  be 
assigned  to  YV  AL. 

Listings  of  example  TF  and  AUXALGO  subroutines  are  given  in  Appendix  A. 
The  examples  given  are  discussed  in  Section  4,  particularly  subsection  4.2.  The 
example  subroutines  use  mnemonic  local  variable  names  to  code  the  algorithms.  This 
facilitates  understanding  and  identification  of  the  algorithm  and  is  a  practice  we 
recommend. 

Once  the  user  has  developed  an  executable  version  of  BETAFACT  which  is  specific 
to  the  algorithm(s)  in  question,  only  one  step  remains  before  interactive  execution 
of  the  problem  can  commence.  This  step  entails  simply  the  collection  of  all  the  pa¬ 
rameter  (both  uncertain  and  independent)  nominal  values  to  be  used  in  the  analyses, 
definition  of  the  K-factor  uncertainties  associated  with  each  uncertain  parameter, 
and  selection  of  the  probability  distribution  type  which  is  to  be  used  in  BETAFACT 
to  model  uncertain  parameter  distributions.  Any  one  of  four  types  of  distributions 
(normal,  lognormal,  Beta,  and  uniform)  may  be  selected.  The  reader  is  cautioned 
that  the  meaning  of  a  K-factor  uncertainty  specification  changes  from  one  probability 
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distribution  type  to  another.  Section  3  should  be  consulted  for  the  precise  meaning 
of  K-factors  as  a  function  of  distribution  type. 

With  the  above  data  and  a  problem  specific  executable  version  of  BETAFACT  in 
hand,  interactive  analysis  of  the  subject  algorithm  may  proceed. 

5.3  INTERACTIVE  EXECUTION. 

In  this  subsection  we  describe  the  specific  meaning  of  prompts  issued  by  BETAFACT 
during  the  course  of  interactive  execution.  The  meaning  of  the  majority  of  these  is 
self-evident  or  has  been  explained,  perhaps  obliquely,  in  earlier  sections  of  this  report. 
We  discuss  them  here  for  completeness  and  to  illustrate  what  is  encountered  during 
the  proper  execution  of  a  BETAFACT  analysis  session. 

We  begin  by  first  mentioning  that  much  of  what  is  entered  interactively  to  BETA- 
FACT  is  not  forever  lost  once  it  scrolls  off  the  user’s  monitor.  Most  input  data 
are  echoed  not  only  to  the  monitor  but  also  to  unit  2  (file  FOR002.DAT  on  VAX 
computers).  This  file  remains  after  execution  of  BETAFACT  is  terminated.  It  can 
be  edited  and  printed  at  will.  Some  data  are  written  to  unit  3  (FOR003.DAT)  only. 
This  data  consists  of  the  x  —  y  pairs,  one  pair  per  line,  which  enable  plotting  of 
computed  and  fitted  algorithm  results  distributions  exterior  to  BETAFACT.  The  file 
also  contains  the  cumulative  probability  distribution  computed  from  the  algorithm 
results  distribution.  The  file  is  generated  only  if  the  user  requests  plot  file  generation 
in  response  to  the  corresponding  program  prompt.  Under  usual  circumstances,  neither 
file  is  excessively  large.  If  BETAFACT  aborts  during  an  interactive  session,  usually 
because  of  an  invalid  input  to  a  prompt,  any  data  previously  written  to  unit  2  or  3 
will  usually  still  be  accessible.  This  is  important  since  no  coding  is  implemented  in 
the  current  version  of  BETAFACT  which  checks  input  as  received  and  enables  the 
user  to  modify  it  if  it  is  incorrect  and  may  lead  to  an  error  termination. 

A  BETAFACT  interactive  session  begins  with  the  user  entering  the  following 
command  (on  VAX  machines)  and  hitting  return: 

RUN  BETAFACT 

The  program  then  begins  to  execute,  issuing  the  following  prompts  (in  some  cases 
we  abbreviate  the  prompts  here)  as  it  proceeds: 

1.  ‘Does  the  algorithm  to  be  evaluated  have  an  overall  uncertainty  associated  with 
its  output?  [Y/NJ’ 

If  it  does,  the  interactive  user  should  type  Y  or  y  and  hit  return.  If  the  basic 
algorithm  does  not  have  associated  uncertainty,  the  user  enters  N  or  n  (or 
anything  else  for  that  matter). 

2.  ‘Are  auxiliary  algorithms  required  to  determine  nominal  parameter  values  which 
are  dependent  on  the  values  of  the  independent  variables?  [Y/N]’ 
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Enter  Y  or  y  if  subroutine  AUXALGO  is  to  be  accessed,  else  enter  N  or  n  (or 
anything  else). 

3.  ‘Enter  the  number  of  parameters  (IP ARAM)  in  the  algorithm  which  have  as¬ 
sociated  uncertainty  (include  in  the  count  whether  or  not  the  algorithm  has  an 
overall  uncertainty),  and  select  the  distribution  type  for  modeling  the  uncer¬ 
tainties  (1,2,3  or  4  for  normal,  lognormal,  Beta,  and  uniform,  respectively).’ 

Enter  the  number  of  parameters  with  associated  uncertainty  (number  of  actual 
parameters  plus  1  if  algorithm  uncertainty  is  to  be  applied)  followed  by  a  space 
and  then  1,2,3,  or  4  to  select  a  distribution  type. 

4.  ‘Enter  the  nominal  values  of  the  parameters  which  are  normally  (lognormally, 
Beta,  or  uniformly)  distributed  and  the  corresponding  K-factors.’ 

On  a  new  line  for  each  uncertain  parameter  enter  its  nominal  value  and  re¬ 
factors  (one  or  two  required  depending  on  distribution  type)  with  each  entry  on 
the  line  separated  by  a  space.  If  two  K-factors  are  entered,  enter  the  low-side 
K-factor  first  and  then  the  high-side  K-factor.  Hit  RETURN  after  each  line  of 
data. 

After  the  entire  set  of  uncertain  parameter  nominal  values  are  entered,  the  data 
are  echoed  to  the  monitor  and  execution  proceeds. 

5.  ‘Enter  the  number  of  independent  variables  ( IXVAL )  which  appear  in  the 
algorithm  or  are  used  in  AUXALGO.' 

Enter  the  number  of  independent  parameters.  At  least  one  must  be  used.  If 
such  parameters  do  not  occur  naturally  in  an  algorithm,  one  may  be  included 
as  an  algorithm  multiplier  with  a  value  of  unity. 

6.  ‘Enter  the  values  of  the  independent  variables.’ 

Enter  the  values,  one  after  another,  separated  by  a  space.  Continue  on  addi¬ 
tional  lines  if  necessary. 

After  the  entire  set  of  independent  variables  are  entered,  the  data  are  echoed 
to  the  monitor. 

7.  ‘Enter  the  number  of  evaluations  of  the  algorithm  to  be  made  for  generating 
the  results  distribution  and  a  negative  seed  integer  for  the  random  number 
generator. 

Enter  a  positive  integer  (say  5000),  a  space,  and  a  negative  integer  (e.g.,  -5731) 
and  then  RETURN.  The  random  number  generator  uses  a  negative  seed  integer. 
If  the  entered  value  is  positive,  the  random  number  generator  uses  its  negative 
value. 

If  parameter  uncertainties  are  to  be  modelled  with  Beta  cumulative  probability 
distributions,  the  following  series  of  prompts  are  then  issued: 

8.  (a)  ‘Select  whether  to  use  tabulated  or  calculated  Beta  cumulative  probability 

distribution  values;  1  -  use  tabulated,  2  -  use  calculated.’ 
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Enter  1  or  2,  as  appropriate. 

If  1  is  entered,  the  program  determines  and  reports  the  number  of  the 
tabulated  distribution  it  will  use  to  model  the  distribution  associated  with 
each  uncertain  parameter. 

If  2  is  entered,  the  program  computes  and  reports  allowable  standard  de¬ 
viation  bound aries  (see  Equation  21)  for  each  uncertain  parameter  and 
prompts. 

(b)  ‘Enter  the  standard  deviations  of  the  Beta  distributed  parameters.’ 

Enter  the  standard  deviation  for  each  parameter,  separating  by  a  space 
and  using  as  many  lines  as  required. 

The  code  echoes  the  input  values  and  reports  the  Beta  cumulative  probabil¬ 
ity  distribution  values  (21  points)  computed  for  each  uncertain  parameter. 

T1  f  code  then  proceeds  to  evaluate  the  algorithm  the  requested  number  of  times,  com¬ 
putes  the  statistics  of  the  results  distribution,  and  then  reports  the  mean,  variance, 
standard  deviation  and  K-factor(s)  determined  for  the  algorithm.  It  then  prompts: 

9.  ‘Do  you  want  plot  files  of  the  results  distribution  and  fit  to  be  generated?  [Y /N].’ 

Enter  Y  or  y  if  x-y  plot  data  are  to  be  saved,  else  enter  N,  n,  or  any  other 
character. 

Finally  the  code  prompts: 

10.  ‘Select  option  to  quit  or  to  run  a  modified  version  of  the  current  problem 
(0,1, 2,..., 16).’ 

A  total  of  17  options  are  available,  one  which  terminates  current  execution 
and  the  other  16  which  enable  various  degrees  of  problem  modification  and  re¬ 
analysis.  See  the  listing  (subroutine  OPTIONS)  for  the  list  of  options  available. 
In  response  to  this  prompt,  type  the  number  of  the  option  selection  and  hit 
RETURN. 

At  this  point,  the  program  reissues  one  or  more  of  the  prompts  described  above  and 
execution  proceeds  as  before.  Thus  the  above  list  of  interactions  of  the  user  with 
BETAFACT  is  complete. 
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SECTION  6 

C  ONCLUSIONS/REC  OMMEND  ATIONS 


In  this  report  we  have  described  a  numerical  approach  for  statistically  combining 
many  sources  of  algorithm  uncertainty  in  order  to  quantify  the  overall  uncertainty  as¬ 
sociated  with  an  algorithm  result.  The  interactive  BETAFACT  program  was  written 
to  accomplish  this  task.  The  code  allows  the  user  to  model  uncertainty  distributions 
using  either  normal,  lognormal,  Beta,  or  uniform  probability  distributions.  We  de¬ 
scribed  in  detail  how  the  code  can  be  used  to  analyze  essentially  any  algorithm  with 
uncorrelated  parameter  uncertainties.  The  only  special  skill  required  by  the  user  to 
use  BETAFACT  to  achieve  this  end  is  moderate  competency  in  writing  FORTRAN 
coding  for  the  algorithms.  During  execution,  the  code  prompts  for  all  data  required 
to  process  the  coded  algorithms. 

We  have  used  BETAFACT  extensively  in  our  LTH-3  Program  lethality  assessment 
activities.  A  code  of  this  type  seems  to  be  essential  in  order  to  prepare  algorithms 
and  define  associated  uncertainties  which  are  then  used  as  input  to  the  DNA  FAST 
code.  BETAFACT  has  also  been  shown  to  be  a  tool  useful  in  program  management 
and  decision  making  processes. 

The  current  version  of  BETAFACT  ha s  some  limitations  which  in  some  appli¬ 
cations  could  be  significant.  We  list  below  some  of  these  limitations  and  gi  ’  our 
recommendations  concerning  whether  or  not  they  should  be  alleviated. 

1 .  The  chief  limitation  of  the  code  is  that  sources  of  uncertainty  are  assumed  to  be 
uncorrelated.  A  substantial  effort  would  be  required  to  implement  a  correlated 
uncertainty  capability  in  the  code.  Such  an  effort  is  recommended  if  there  is  a 
desire  to  enhance  the  generality  of  the  code. 

2.  A  possibly  important  limitation  is  the  requirement  that  all  sources  of  uncer¬ 
tainty  in  a  given  analysis  must  be  modelled  with  one  type  of  probability  dis¬ 
tribution  (e.g.,  all  normal  or  all  lognormal,  etc.).  It  would  be  straightforward 
to  modify  BETAFACT  to  allow  modelling  different  sources  of  uncertainty  with 
different  probability  distributions  in  the  same  analysis.  The  resulting  algorithm 
results  distributions  would  tend  to  be  hybrid  in  this  case  and  perhaps  not  easily 
or  accurately  modelled  by  one  of  the  standard  types  of  distributions.  If  a  tab¬ 
ular  form  of  the  algorithm  results  distribution  is  useful  (for  instance,  if  FAST 
was  modified  to  accept  and  use  arbitrary  probability  distributions),  then  mod¬ 
ification  of  BETAFACT  to  allow  mixing  of  distribution  types  in  an  analysis  is 
recommended. 

3.  Only  four  types  of  probability  distributions  (normal,  lognormal,  Beta,  and  gen¬ 
eralized  uniform)  are  available  in  BETAFACt.  We  recommend  adding  addi¬ 
tional  distribution  types  to  the  code.  The  effort  involved  would  not  be  substan¬ 
tial. 
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APPENDIX  A 

EXAMPLE  USER  SUBROUTINES 


1.  This  is  the  simple  TF  subroutine  used  in  the  verification  analyses  of  the  BETAFACT 
code.  For  a  discussion  of  the  simple  algorithm  and  the  analysis  results,  see  subsection 
4.1  of  the  BETAFACT  documentation. 


SUBROUTINE  TF ( ANOM , A , XVAL , YVAL) 

C 

C  THIS  EXAMPLE  FOR  SUBROUTINE  TF  CONTAINS  THE  FORTRAN  CODING  FOR 

C  THE  SIMPLE  ALGORITHM  CONSIDERED  IN  SUBSECTION  4.1  OF  THE  BETAFACT 

C  CODE  DOCUMENTATION : 

C 

C  R  =  X 

C 

C  THE  VARIABLE  X,  WHICH  HAS  ASSOCIATED  UNCERTAINTY,  HAS  A 

C  RANDOMIZED  VALUE  WHICH  IS  PASSED  TO  THIS  SUBROUTINE  AS  THE 

C  FIRST  ELEMENT,  A(l),  IN  ARRAY  A.  THE  LOCAL  NAME  FOR  THIS 

C  VARIABLE  IS  ADUM.  THE  RANDOMIZED  VARIABLE  VALUE  IS  MULTIPLIED 

C  BY  THE  VALUE  OF  THE  INDEPENDENT  VARIABLE  (WHICH  HAS  VALUE  1.0) 

C  WHICH  IS  STORED  AS  THE  FIRST  ELEMENT  OF  ARRAY  XVAL.  THE 

C  OUTPUT  R  OF  THE  SIMPLE  ALGORITHM  IS  ASSIGNED  TO  YVAL,  AS 

C  REQUIRED  BY  BETAFACT. 

C 

REAL*4  ANOM(*) ,  A(*),  XVAL(*) 

ADUM  =  A ( 1 ) 

YVAL  =  XVAL(l) *ADUM 

RETURN 

END 
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2.  The  following  are  the  AUXALGO  and  TF  subroutines  used  in  the  realistic 
algorithm  example  analyses  described  in  subsection  4.2  of  the  BETAFACT  documen¬ 
tation. 


SUBROUTINE  AUXALGO (ANOM ,  XVAL) 

C 

C  THIS  SAMPLE  SUBROUTINE  CONTAINS  CODING  NEEDED  FOR  EVALUATING  THE 

C  MODIFIED  VERSION  OF  THE  HKL  ALGORITHM  DESCRIBED  IN  SUBSECTION  4.2 

C  OF  THE  BETAFACT  DOCUMENTATION. 

C 

REAL*4  AN0M(*) ,  XVAL(*) 

AI  =  XVAL(l) 

AH  =  XVALC2) 

IF  (AI.LT.3.5)  THEN 
AN0M(3)  =  AI 
AN0M(4)  =  AH 
RETURN 
ELSE 

FAC  =  SqRT(2.0*AH) 

AN0M(3)  =  (3.5  +  0.08*(AI  -  3.5))*FAC 
ANOM(4)  =  AH  -  (0.05  +  0.001*(AI  -  3.5))*FAC 
RETURN 
END  IF 
RETURN 
END 

SUBROUTINE  TF (ANOM , A , XVAL ,  YVAL) 

C 

C  THIS  SAMPLE  SUBROUTINE  CONTAINS  CODING  FOR  THE  MODIFIED  VERSION 

C  OF  THE  HKL  ALGORITHM  DESCRIBED  IN  SUBSECTION  4.2  OF  THE  BETAFACT 

C  DOCUMENTATION . 

C 

REAL*4  AN0M(*) ,  A(*),  XVAL(*) 

ARHO  =  A(l) 

ASIG  =  A(2) 

AIR  =  A(3) 

AHR  *  A (4) 

YVAL  =  0.5*AIR*AIR/(ARH0*ASIG*AHR**1.5) 

YVAL  =  YVAL*10**6 . 

RETURN 

END 
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APPENDIX  B 

BETAFACT  CODE  LISTING 


PROGRAM  BETAFACT 
C 

C  PROGRAM  BETAFACT  IS  USED  TO  COMPUTE  THE  COMBINED  UNCERTAINTY 

C  TO  BE  ASSOCIATED  WITH  AN  ALGORITHM  WHICH  IS  A  FUNCTION  OF 

C  SEVERAL  PARAMETERS,  EACH  OF  WHICH  MAY  HAVE  ASSOCIATED 
C  UNCERTAINTY.  THE  ALGORITHM,  ITSELF,  MAY  HAVE  SPECIFIED 

C  UNCERTAINTY  ASSOCIATED  WITH  ITS  OUTPUT,  REGARDLESS  OF 

WHETHER  OR  NOT  ITS  PARAMETERS  HAVE  UNCERTAIN  VALUES.  ALL 
THESE  SOURCES  OF  UNCERTAINTY  ARE  ACCOUNTED  FOR  IN  BETAFACT 
C  IN  ORDER  TO  OBTAIN  THE  FINAL  UNCERTAINTY  ESTIMATE  FOR  THE 

C  EVALUATED  ALGORITHM.  ALL  UNCERTAINTIES  USED  IN  THE  PROGRAM 

C  ARE  SPECIFIED  IN  THE  FORM  OF  EQUIVALENT  K-FACTORS. 

C 

C  THE  UNCERTAINTIES  ASSOCIATED  WITH  PARAMETERS  AND  ALGORITHMS 

C  MAY  BE  MODELLED  AS  NORMALLY,  LQGNORMALLY,  BETA,  OR  UNIFORMLY 

C  DISTRIBUTED  IN  BETAFACT.  THE  USER  IS  PROMPTED  TO  DEFINE, 

C  INTERACTIVELY,  THE  REQUIRED  DISTRIBUTION  TYPES  AND,  IN  FACT, 

C  THE  ENTIRE  SPECIFICATION  OF  THE  PROBLEM  TO  BE  SOLVED.  IT  IS 
C  ONLY  REQUIRED  THAT  THE  USER  PROVIDE  A  SUBROUTINE  (CALLED  TF) 

C  TO  THE  PROGRAM  WHICH  CONTAINS  THE  CODING  FOR  THE  ALGORITHM 

C  WHICH  IS  TO  BE  EVALUATED.  THE  USER  MAY  ALSO  PROVIDE  AN 

C  AUXILIARY  SUBROUTINE  (CALLED  AUXALGO)  TO  EVALUATE  THE  VALUES 

C  OF  PARAMETERS  USED  IN  THE  MAIN  ALGORITHM  WHICH  ARE  OBTAINABLE 

C  FORM  OTHER  SIMPLE  ALGORITHMS  OR  RELATIONS. 

C 

C 

C  BETAFACT  WAS  DEVELOPED  ON  A  VAX  TYPE  COMPUTER  OPERATING  UNDER  VMS. 

C  TO  USE  BETAFACT,  THE  USER  FIRST  MUST  COMPILE  THE  SUBROUTINE  TF 

C  (AND  SUBROUTINE  AUXALGO  IF  THERE  IS  ONE)  AND  THEN  LINK  BOTH  TO 

C  THE  BETAFACT  OBJECT  FILE  AS  FOLLOWS: 

C 

C  LINK  BET AF ACT, TF, AUXALGO (if  there  is  one) 

C 

C  THE  RESULT  OF  THIS  PROCESS  IS  A  FILE  CALLED  BETAFACT.EXE,  WHICH 

C  IS  THE  BETAFACT  EXECUTABLE  FILE.  ONCE  THIS  FILE  IS  CREATED, 

C  ITS  EXECUTION  IS  ACHIEVED  SIMPLY  BY  TYPING 

C 

C  RUN  BETAFACT 

C 

C  HITTING  RETURN,  AND  THEN  ANSWERING  THE  PROMPTS  ISSUED  BY  THE 

C  PROGRAM . 
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A  SIGNIFICANT  PORTION  OF  BETAFACT  IS  FLEXIBLY  DIMENSIONED. 

ONLY  THE  ARRAYS  EXPLICITLY  RELATED  TO  THE  LIBRARY  OF  81 
BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS  USED  BY  THE 
CODE  HAVE  FIXED  DIMENSIONS.  IF  THE  NUMBER  OF  LIBRARY 
BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS  NEEDS  TO  BE 
INCREASED,  THE  DIMENSIONS  OF  THE  FOLLOWING  ARRAYS  WILL 
HAVE  TO  BE  CHANGED:  BD,  DBA,  DBB,  AMED,  ALPH,  AND  BET. 

ALL  OTHER  ARRAYS  USED  BY  THE  PROGRAM  ARE  DIMENSIONED  VIA 
A  PARAMETER  STATEMENT  WHICH  DEFINES  THE  FOLLOWING  INTEGER 
VARIABLES : 

MAXPRM  -  MAXIMUM  NUMBER  OF  VARIABLES  WITH  UNCERTAINTY. 

MAXVAR  -  MAXIMUM  NUMBER  OF  INDEPENDENT  PARAMETERS. 

MAXSIZ  -  MAXIMUM  NUMBER  OF  ALGORITHM  EVALUATIONS. 

MAXINT  -  MAXIMUM  NUMBER  OF  INTERVALS  USED  TO  CALCULATE  BETA 
CUMULATIVE  PROBABILITY  DISTRIBUTIONS. 

MAXIN1  -  MAXIMUM  NUMBER  OF  INTERVALS  USED  TO  DEFINE  HISTOGRAMS. 

OTHER  PRINCIPAL  INTEGER  VARIABLES  AND  CONTROL  PARAMETERS 
USED  BY  THE  CODE  ARE  AS  FOLLOWS: 

IPARAM  -  NUMBER  OF  PARAMETERS  WITH  UNCERTAINTY. 

IPARA  -  EQUAL  TO  IPARAM  IF  ALGORITHM  DOES  NOT  HAVE  AN  OVERALL 
SPECIFIED  UNCERTAINTY;  ELSE  IPARA  =  IPARAM  -  1. 

ITYPE  -  DISTRIBUTION  TYPE  USED  TO  MODEL  UNCERTAINTIES: 

=  1  NORMAL  DISTRIBUTION. 

*  2  LOGNORMAL  DISTRIBUTION. 

=  3  BETA  DISTRIBUTION. 

*  4  UNIFORM  DISTRIBUTION. 

ISIZE  -  NUMBER  OF  ALGORITHM  EVALUATIONS  TO  BE  MADE. 

I SEED  -  SEED  (NEGATIVE  AND  ODD)  USED  TO  INITIALIZE  THE  RANDOM 
NUMBER  GENERATOR. 

I ALGO  -  FLAG  INDICATING  SPECIFIED  ALGORITHMIC  UNCERTAINTY: 

*  1  ALGORITHM  HAS  OVERALL  UNCERTAINTY. 

=  ANYTHING  ELSE;  DOESN’T  HAVE  OVERALL  UNCERTAINTY. 
IXVAL  -  NUMBER  OF  INDEPENDENT  PARAMETERS  USED  IN  ALGORITHM. 

THESE  DO  NOT  HAVE  ASSOCIATED  UNCERTAINTY. 

ITABL  -  FLAG  DEFINING  WHETHER  TABULATED  OR  CALCULATED  BETA 
CUMULATIVE  PROBABILITY  DISTRIBUTIONS  ARE  TO  BE  USED: 

=  1  USE  TABULATED  DISTRIBUTIONS. 

=  2  USE  CALCULATED  DISTRIBUTIONS. 

IXALG  -  FLAG  WHICH  SPECIFIES  WHETHER  OR  NOT  AUXILIARY 
ALGORITHMS  IN  USER  SUPPLIED  SUBROUTINE  AUXALGO 
ARE  TO  BE  ACCESSED:  IF  »  1,  THEN  SUBROUTINE  IS 
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ACCESSED,  OTHERWISE  NOT  ACCESSED. 

IOPT  -  ANALYSIS  MODIFICATION  OPTIONS  FLAG.  SEE  SUBROUTINE 
OPTIONS  FOR  AVAILABLE  OPTIONS. 

NINTV  -  NUMBER  OF  INTERVALS  USED  TO  CALCULATE  BETA 

CUMULATIVE  PROBABILITY  DISTRIBUTIONS.  RECOMMENDED 
NUMBER  IS  DEFAULT  NINTV  =  MAXINT. 

NINC  -  NUMBER  OF  INTERVALS  USED  TO  GENERATE  HISTOGRAMS. 
RECOMMENDED  NUMBER  IS  DEFAULT  NINC  =  MAXIN1 . 

NINP  -  UNIT  NUMBER  FOR  INTERACTIVE  INPUT  (DEFAULT  =  5) . 

NOUT  -  UNIT  NUMBER  FOR  TERMINAL  SCREEN  OUTPUT  (DEFAULT  =  6) . 

NPRT  -  UNIT  NUMBER  FOR  PRINTABLE  OUTPUT  (DEFAULT  =  2) . 

NPLT  -  UNIT  NUMBER  FOR  PLOT  FILE  (HISTOGRAM)  OUTPUT  (DEFAULT  =  3) . 

THE  PRINCIPAL  REAL  ARRAYS  AND  VARIABLES  USED  BY  THE  PROGRAM  ARE: 

ANOM  -  ARRAY  OF  NOMINAL  VALUES  OF  PARAMETERS  WITH  UNCERTAINTY. 

AKLO  -  ARRAY  OF  PARAMETER  LOW-SIDE  K-FACTOR  UNCERTAINTIES. 

AKHI  -  ARRAY  OF  PARAMETER  HIGH-SIDE  K-FACTOR  UNCERTAINTIES. 

ALO  -  ARRAY  OF  PARAMETER  LOW-SIDE  EXTREME  VALUES. 

AHI  -  ARRAY  OF  PARAMETER  HIGH-SIDE  EXTREME  VALUES. 

A  -  ARRAY  OF  RANDOM  VALUES  OF  PARAMETERS  FOR  A  SINGLE 

ALGORITHM  EVALUATION  LOOP. 

XVAL  -  ARRAY  OF  INDEPENDENT  PARAMETER  VALUES. 

YVAL  -  RESULT  OBTAINED  IN  SINGLE  EVALUATION  OF  ALGORITHM. 

RESULT  -  ARRAY  OF  ALL  COMPUTED  YVAL  RESULTS. 

ASIG  -  ARRAY  OF  SPECIFIED  BETA  DISTRIBUTION  STANDARD  DEVIATIONS; 
USED  TO  CALCULATE  BETA  DISTRIBUTIONS. 

BD  -  ARRAY  OF  LIBRARY  OR  CALCULATED  BETA  CUMULATIVE 
PROBABILITY  DISTRIBUTIONS. 

ALPHA  -  BETA  DISTRIBUTION  ALPHA  PARAMETER. 

BETA  -  BETA  DISTRIBUTION  BETA  PARAMETER. 

ALPH  -  ARRAY  OF  ALPHA  PARAMETER  VALUES  USED  TO  GENERATE  THE 
LIBRARY  OF  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS. 

BET  -  ARRAY  OF  BETA  PARAMETER  VALUES  USED  TO  GENERATE  THE 
LIBRARY  OF  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS. 

AMED  -  ARRAY  OF  MEDIAN  (50*/.)  VALUES  OF  THE  LIBRARY  OF 
BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS. 

RN  -  VARIABLE  EQUAL  TO  THE  MOST  RECENT  RANDOM  NUMBER. 

AMEAN  -  MEAN  OF  THE  ALGORITHM  RESULTS  DISTRIBUTION. 

VAR  -  VARIANCE  OF  THE  ALGORITHM  RESULTS  DISTRIBUTION. 

SIGMA  -  STANDARD  DEVIATION  OF  ALGORITHM  RESULTS  DISTRIBUTION. 

AKEST  -  COMBINED  EFFECTIVE  K-FACTOR  ESTIMATE  FOR  CASES  OF 
NORMAL  AND  LOGNORMAL  MODELLED  UNCERTAINTIES. 

AKEST 1  -  LOW-SIDE  COMBINED  EFFECTIVE  K-FACTOR  ESTIMATE  FOR 
BETA  AND  UNIFORM  MODELLED  UNCERTAINTIES. 

AKEST2  -  HIGH-SIDE  COMBINED  EFFECTIVE  K-FACTOR  ESTIMATE  FOR 
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BETA  AND  UNIFORM  MODELLED  UNCERTAINTIES. 

AKAVE  -  K-FACTOR  EqUAL  TO  THE  AVERAGE  OF  AKEST1  AND  AKEST2 . 

THE  PRINCIPAL  INTEGER  ARRAYS  USED  BY  THE  PROGRAM  ARE: 

IBD  -  ARRAY  OF  THE  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTION 
NUMBER  ASSOCIATED  WITH  EACH  BETA  DISTRIBUTED  PARAMETER. 
IFREQ  -  FREQUENCY  OF  OCCURRENCE  ARRAY  USED  TO  ACCUMULATE 
HISTOGRAM  DATA. 

PARAMETER  (MAXPRM=20,  MAXVAR*10,  MAXSIZ* lOOOO, 

1  MAXINT*100 ,  MAXIN1M00) 

COMMON  /BLK01/  IPARAM,  I TYPE,  ISIZE,  ISEED,  IALGO ,  IXVAL,  ITABL, 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

COMMON  /BLK02/  BD(21,81) 

REALM  ANOM(MAXPRM) ,  ALO(MAXPRM) ,  AHI (MAXPRM) 

REALM  AKLO (MAXPRM),  AKHI (MAXPRM) ,  A (MAXPRM) 

REALM  XVAL(MAXVAR) ,  RESULT(MAXSIZ) 

REALM  ASIG(2, MAXPRM),  AREA (MAXINT),  Y(MAXINT+1) 

REALM  XINT(MAXIN1+1)  ,  XOUT(MAXINl) 

REALM  XXLN(MAXIN1+1)  ,  XXII(MAXINl) 

INTEGERM  IBD  (MAXPRM)  ,  IFREQ  (MAXIN1) 

CONTROL  PARAMETER  'IOPT’  INITIALLY  SET  EQUAL  TO  16 
SINCE  A  COMPLETELY  NEW  PROBLEM  IS  TO  BE  DEFINED. 


C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 


IOPT  ■  16 
NINC  =  MAXINT 
NINTV  *  MAXIN1 

ASSIGN  DEFAULT  UNIT  NUMBERS  FOR  INPUT,  OUTPUT,  ETC: 

NINP  -  UNIT  NUMBER  FOR  INPUT  FROM  TERMINAL. 
NOUT  -  UNIT  NUMBER  FOR  OUTPUT  TO  TERMINAL. 
NPLT  -  UNIT  NUMBER  FOR  PLOT  DATA. 

NPRT  -  UNIT  NUMBER  FOR  HARDCOPY  OUTPUT. 

NINP  =  5 
NOUT  =»  6 
NPLT  *  3 
NPRT  =  2 

ASSIGN  DEFAULT  VALUES  FOR  OTHER  CONTROL  PARAMETERS. 

IALGO  =  0 
IXALG  =  0 
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IXVAL  *  1 
ITABL  *  1 
C 

C  CALL  SUBROUTINE  CONTRL  TO  INPUT  PROBLEM  DEFINITION  AND 
MAIN  CONTROL  PARAMETERS. 

10  CALL  CNTRL(ANOM,  ALO,  AHI,  AKLO,  AKHI,  XVAL,  ASIG,  Y,  AREA,  IBD) 

CALL  SUBROUTINE  SOLVE  TO  EVALUATE  THE  PROBLEM  ALGORITHM. 

CALL  SOLVE(ANOM,  ALO,  AHI,  AKLO,  AKHI,  XVAL,  A,  IBD,  RESULT) 

CALL  SUBROUTINE  SORT  TO  ORDER  THE  ALGORITHM  RESULTS 
C  DISTRIBUTION  FROM  LOWEST  TO  HIGHEST  VALUE. 

C 

CALL  SORT(RESULT,ISIZE) 

C 

C  CALL  SUBROUTINE  STATISTIC  TO  EVALUATE  STATISTICS  FOR  THE 

C  ORDERED  RESULTS  DISTRIBUTION,  TO  ESTIMATE  OVERALL  ALGORITHM 

C  UNCERTAINTY,  AND  TO  GENERATE  PLOT  FILES  CONTAINING  HISTOGRAM 

C  DATA  AND  FITS  TO  THE  HISTOGRAM  DATA. 

C 

CALL  STATISTIC (RESULT,  XINT,  XOUT,  XXLN,  XXII,  I FREQ) 

C 

C  CALL  SUBROUTINE  OPTIONS  TO  MODIFY  THE  CURRENT  PROBLEM 

C  DEFINITION  WITHOUT  EXITING  FROM  BETAFACT. 

C 

CALL  OPTIONS 
C 

C  TERMINATE  THE  PRESENT  ANALYSIS  ONLY  IF  IOPT  =  0. 

C 

IF  (IOPT.NE.O)  GO  TO  10 
END 
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SUBROUTINE  CNTRL(ANOM,  ALO,  AHI,  AKLO,  AKHI,  XVAL,  ASId,  Y, 

1  AREA,  IBD) 

SUBROUTINE  CNTRL  CONTROLS  THE  INPUT  OF  THE  DATA  WHICH  DEFINE  A 
PROBLEM  OR  MODIFICATIONS  TO  A  PREVIOUSLY  DEFINED  PROBLEM.  THE 
SUBROUTINE  INTERACTIVELY  PROMPTS  THE  USER  TO  SPECIFY  THE  VALUES 
OF  THE  FOLLOWING  CONTROL  PARAMETERS  AND  INTEGER  VARIABLES: 

I ALGO  -  FLAG  IDENTIFYING  USE  OF  OVERALL  ALGORITHM  UNCERTAINTY. 
IXALG  -  FLAG  CONTROLLING  ACCESSING  OF  AUXILIARY  ALGORITHMS. 
IPARAM  -  NUMBER  OF  PARAMETERS  WITH  UNCERTAINTY. 

I TYPE  -  DISTRIBUTION  TYPE  FOR  MODELLING  UNCERTAINTIES. 

ISIZE  -  NUMBER  OF  REQUESTED  ALGORITHM  EVALUATIONS. 

ISEED  -  SEED  INTEGER  (<0  AND  ODD)  FOR  RANDOM  NUMBER  GENERATOR. 

WHETHER  OR  NOT  ONE  OR  ALL  OF  THE  ABOVE  INTEGER  VARIABLES  NEEDS 
TO  BE  ENTERED  IS  DETERMINED  BY  THE  CURRENT  VALUE  OF  THE  OPTIONS 
CONTROL  PARAMETER  I OPT. 

SUBROUTINE  CNTRL  ALSO  CALLS  TWO  SUBROUTINES: 

SUBROUTINE  INPUT  -  FOR  ENTERING  PARAMETER  NOMINAL  VALUES, 

ASSOCIATED  UNCERTAINTIES,  AND  THE  VALUES 
OF  INDEPENDENT  PARAMETERS  WHICH  DO  NOT 
HAVE  ASSOCIATED  UNCERTAINTY. 

SUBROUTINE  ASSIGNBETA  -  WHICH  CONTROLS  ASSOCIATION  OF  A  SPECIFIC 

TABULATED  OR  COMPUTED  BETA  CUMULATIVE 
PROBABILITY  DISTRIBUTION  WITH  PARAMETERS 
HAVING  BETA  DISTRIBUTED  UNCERTAINTY. 

COMMON  /BLKOl/  IPARAM,  I TYPE,  ISIZE,  ISEED,  I ALGO ,  IXVAL,  ITABL, 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

REAL*4  ANOMO) ,  AL0(*),  AHI(*),  AKLQ(*)  ,  AKHI(*),  XVAL(*) 

REAL*4  ASIG(2,*) ,  Y(*) ,  AREA(*) 

INTEGERS  IBD(*) 

CHARACTER* 1  ANSI 

THE  VALUE  OF  THE  OPTION  CONTROL  PARAMETER  'IOPT’ 

EQUALS  16  ONLY  IF  THIS  IS  A  NEW  PROBLEM  OR  THE 
CURRENT  PROBLEM  IS  TO  BE  COMPLETELY  MODIFIED. 

FOR  IOPT  *  16,  THE  CONTROL  PARAMETERS  ‘I ALGO'  AND 
'IXALG'  NEED  TO  BE  SPECIFIED. 

IF  (IOPT. NE. 16)  GO  TO  10 

SPECIFY  CONTROL  PARAMETER  'IALGO'  WHICH  FLAGS  WHETHER  OR  NOT 
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C  THE  ALGORITHM  HAS  AN  OVERALL  ASSOCIATED  UNCERTAINTY. 

C  IF  I ALGO  *  1,  THEN  THE  ALGORITHM  HAS  OVERALL  UNCERTAINTY. 

C  OTHERWISE  IT  DOES  NOT. 

C 

WRITE (NOUT, 1000) 

1000  FORMAT (/, '  DOES  THE  ALGORITHM  TO  BE  EVALUATED  HAVE  AN  OVERALL '  ,  / , 

1  '  UNCERTAINTY  ASSOCIATED  WITH  ITS  OUTPUT?  [Y/N]',/) 

READ(NINP, 1001)  ANSI 

1001  FORMAT(Al) 

IF  ((ANS1.EQ.  'Y') .OR. (ANS1.EQ.  'y'))  I ALGO  =  1 


SPECIFY  CONTROL  PARAMETER  ‘IXALG'  WHICH  FLAGS  WHETHER  OR  NOT 
C  THE  ALGORITHM  MUST  ACCESS  AUXILIARY  ALGORITHMS  TO  OBTAIN  VALUES 

C  FOR  SOME  OF  ITS  PARAMETERS.  IF  IXALG  =  1,  THEN  AUXILIARY 

C  ALGORITHMS,  CODED  IN  THE  USER  SUPPLIED  SUBROUTINE  AUXALGO, 

C  WILL  BE  ACCESSED.  OTHERWISE  AUXILIARY  ALGORITHMS  ARE  NOT 

C  USED . 

C 

WRITE(NOUT.llOO) 

1100  FORMAT (/, ’  ARE  AUXILIARY  ALGORITHMS  REQUIRED  TO  DETERMINE',/, 

1  '  NOMINAL  PARAMETER  VALUES  WHICH  ARE  DEPENDENT  ON',/, 

2  '  THE  VALUES  OF  THE  INDEPENDENT  VARIABLES?  [Y/N]',/) 
READ (NINP, 1001)  ANSI 

IF  ((ANSI .EQ.  ' Y ' ) . OR . (ANS 1 . EQ .  'y'))  IXALG  *  1 
C 

C  IF  THIS  IS  A  MODIFICATION  ANALYSIS,  CHECK  THE  VALUE  OF  'IOPT' 

C  TO  DETERMINE  IF  THE  NUMBER  OF  ALGORITHM  PARAMETERS  AND/OR 

C  DISTRIBUTION  TYPE  SPECIFICATION  NEED  TO  BE  CHANGED. 

C 

10  IF  ( (IOPT . Eq . 1 ) . OR. (IOPT . Eq . 2) . OR . (IOPT . EQ . 3) . OR . (IOPT . EQ . 5) . OR . 
1  (IOPT. EQ. 6). OR. (IOPT. EQ. 8). OR. (IOPT. EQ. 11))  GO  TO  20 
C 


C  SPECIFY  THE  NUMBER  ‘IPARAM'  OF  ALGORITHM  PARAMETERS  WHICH 

C  HAVE  ASSOCIATED  UNCERTAINTY  AND  THE  DISTRIBUTION  TYPE  ' ITYPE ’ 
C  TO  BE  USED  TO  MODEL  THE  PARAMETER  UNCERTAINTIES. 

C 


WRITE (NOUT, 1200) 

1200  FORMAT (/ , '  ENTER  THE  NUMBER  OF  PARAMETERS  (IPARAM)  IN  THE',/, 

1  '  ALGORITHM  WHICH  HAVE  ASSOCIATED  UNCERTAINTY  ',/, 

2  '  (INCLUDE  IN  THE  COUNT  WHETHER  OR  NOT  THE  ALGORITHM 

3  '  HAS  AN  OVERALL  UNCERTAINTY),  AND  SELECT  THE  ',/, 

4  '  DISTRIBUTION  TYPE  FOR  MODELLING  THE  UNCERTAINTIES  ’ 

5  '  FROM  THE  FOLLOWING  (ENTER  1,  2,  3,  OR  4):',/,/, 

6  '  1  -  NORMAL  DISTRIBUTION.',/, 

7  '  2  -  LOGNORMAL  DISTRIBUTION.',/, 

8  ’  3  -  BETA  DISTRIBUTION.',/, 


,/, 


,/, 
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9  ’  4  -  UNIFORM  DISTRIBUTION. 1 ,/) 

READ(NINP,*)  IPARAM,  I TYPE 
20  IF  (I0PT.EQ.2)  GO  TO  30 
C 

C  CALL  SUBROUTINE  INPUT  TO  ENTER  DEPENDENT  PARAMETER  NOMINAL  VALUES, 
C  ASSOCIATED  K-FACTOR  UNCERTAINTIES,  AND  THE  VALUE(S)  OF  THE 

C  INDEPENDENT  PARAMETER(S)  FOR  WHICH  THE  ALGORITHM  IS  TO  BE 

C  EVALUATED. 

C 

CALL  INPUT (ANOM,  ALO,  AHI ,  AKLO,  AKHI,  XVAL) 

C 

C  IF  THIS  IS  A  MODIFICATION  ANALYSIS,  CHECK  THE  VALUE  OF  'IOPT' 

C  TO  DETERMINE  IF  THE  NUMBER  'ISIZE1  OF  EVALUATIONS  OF  THE 

C  ALGORITHM  AND/OR  THE  SEED  INTEGER  'ISEED1  FOR  THE  RANDOM 

C  NUMBER  GENERATOR  NEED  TO  BE  CHANGED. 

C 

30  IF  ((IOPT.Eq.l) .OR. (IOPT. Eq. 3) .OR. (IOPT. Eq. 4) .OR. (IOPT. Eq. 6) .OR. 

1  (IOPT. Eq. 7) .OR. (IOPT. Eq. 10) .OR. (IOPT. Eq. 13))  GO  TO  40 
C 

C  SPECIFY  THE  NUMBER  'ISIZE1  OF  ALGORITHM  EVALUATIONS  TO  BE 
C  PERFORMED  AND  A  SEED  INTEGER  'ISEED'  FOR  THE  RANDOM  NUMBER 

C  GENERATOR. 

C 

WRITE (NOUT, 1300) 

1300  FORMAT(/,'  ENTER  THE  NUMBER  OF  EVALUATIONS  (ISIZE)  OF  THE1,/, 

1  1  ALGORITHM  TO  BE  MADE  FOR  GENERATING  THE  RESULTS1,/, 

2  1  DISTRIBUTION  AND  A  NEGATIVE  SEED  INTEGER  (ISEED)1,/, 

3  1  FOR  THE  RANDOM  NUMBER  GENERATOR.1,/) 

READ(NINP,*)  ISIZE,  ISEED 

40  CONTINUE 
C 

C  CALL  SUBROUTINE  ASSGNBETA  IF  THE  PRESENT  ANALYSIS  INVOLVES 

C  MODELLING  UNCERTAINTIES  WITH  CUMULATIVE  BETA  DISTRIBUTIONS. 

C 

IF  ((IOPT.Eq.l) .OR. (IOPT. Eq. 2) .OR. (IOPT. Eq. 5))  GO  TO  50 
IF  (ITYPE.EQ.3)  THEN 

CALL  ASSGNBETA (ANOM,  ALO,  AHI,  AKLO,  AKHI,  ASIG,  Y,  AREA,  IBD) 
ENDIF 

50  CONTINUE 
RETURN 
END 


SUBROUTINE  INPUT (ANOM,  ALO,  AHI,  AKLO,  AKHI,  XVAL) 


C 

C  SUBROUTINE  INPUT  INTERACTIVELY  READS  USER  SUPPLIED  DATA  WHICH 

C  DEFINE  THE  NOMINAL  VALUES  OF  ALGORITHM  PARAMETERS  WITH  UNCERTAINTY, 

C  THE  ASSOCIATED  K-FACTOR  UNCERTAINTIES,  AND  THE  VALUES  OF  INDEPENDENT 

C  PARAMETERS  WHICH  DO  NOT  HAVE  ASSOCIATED  UNCERTAINTY.  THE  NUMBER 

C  OF  THE  LATTER  TYPE  OF  PARAMETERS,  IXVAL,  IS  PROMPTED  FOR  BY  THE 
C  SUBROUTINE.  THE  PARAMETER  NOMINAL  VALUES  AND  K-FACTOR  UNCERTAINTIES 

C  ARE  USED  IN  THE  SUBROUTINE  TO  COMPUTE  THE  EXTREME  (LOW  AND  HIGH) 

C  VALUES  CORRESPONDING  TO  EACH  PARAMETER  WITH  UNCERTAINTY. 

C 

C  THE  PRINCIPAL  ARRAYS  DEFINED  IN  SUBROUTINE  INPUT  ARE  AS  FOLLOWS: 

C 

C  ANOM  -  NOMINAL  VALUES  OF  THE  PARAMETERS  WITH  UNCERTAINTY. 

C  AKLO  -  LOW-SIDE  K-FACTOR  UNCERTAINTIES  FOR  THE  PARAMETERS. 

C  AKHI  -  HIGH-SIDE  K-FACTOR  UNCERTAINTIES  FOR  THE  PARAMETERS. 

C  ALO  -  LOW-SIDE  OF  THE  RANGE  FOR  EACH  PARAMETER. 

C  AHI  -  HIGH-SIDE  OF  THE  RANGE  FOR  EACH  PARAMETER. 

C  XVAL  -  VALUES  OF  THE  INDEPENDENT  PARAMETERS  (NO  UNCERTAINTY) . 

C 

COMMON  /BLKOl/  IPARAM,  I TYPE,  ISIZE,  I SEED,  IALGO ,  IXVAL,  ITABL , 

1  IXALG ,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

REAL*4  ANOM(*),  ALO(*),  AHI(*),  AKLO(*) ,  AKHI(*) ,  XVAL(*) 

C 

C  DETERMINE  IF  IOPT  VALUE  REQUIRES  INPUT  OF  PARAMETER  DATA. 

C 

IF  ((IOPT.EQ. 1) .OR. (IOPT.EQ . 2) .OR. (IOPT. EQ .5))  GO  TO  90 
C 

C  ENTER  DATA  FOR  PARAMETERS  WHICH  ARE  NORMALLY  DISTRIBUTED. 

C 

IF  (ITYPE.EQ. 1)  THEN 
WRITE(NOUT, 1000) 

DO  10  I  *  1,  IPARAM 

READ (NINP,*)  ANOM (I) ,  AKLO(I) 

AHI  (I)  =>  ANQM(I)*AKLO(I) 

ALO (I)  =  ANOM ( I )  -  (AHI (I)  -  ANOM (I)) 

10  CONTINUE 

C 

C  ECHO  PARAMETER  INPUT  DATA  AND  COMPUTED  RANGE. 

C 

WRITE(NPRT.llOO) 

WRITE(NQUT, 1 100) 

DO  20  I  =  1,  IPARAM 

WRITE (NPRT, 1200)  I,  ANOM(I),  AKLO(I),  ALO(I),  AHI (I) 
WRITE(NOUT, 1200)  I,  ANOM(I) ,  AKLO(I) ,  ALO(I),  AHI(I) 
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20  CONTINUE 
C 

C  ENTER  DATA  FOR  PARAMETERS  WHICH  ARE  LOGNORMALLY  DISTRIBUTED. 

C 

ELSE  IF  (ITYPE.EQ.2)  THEN 
WRITE(NOUT, 1300) 

DO  30  I  *  1,  IPARAM 

READ(NINP,*)  ANOM(I),  AKLO(I) 

ALO(I)  *  ANOM(I)/AKLO(I) 

AHI(I,  *  ANOM(I)*AKLO(I) 

30  CONTINUE 

C 

C  ECHO  PARAMETER  INPUT  DATA  AND  COMPUTED  RANGE. 

C 

WRITE(NPRT, 1400) 

WRITE(NOUT, 1400) 

DO  40  I  *  1,  IPARAM 

WRITE(NPRT, 1500)  I,  ANOM(I) ,  AKLO(I),  ALO(I),  AHI(I) 
WRITE(N0UT,15O0)  I,  ANOM(I) ,  AKLO(I),  ALO(I),  AHI(I) 

40  CONTINUE 

C 

C  ENTER  DATA  FOR  PARAMETERS  WHICH  ARE  BETA  DISTRIBUTED. 

C 

ELSE  IF  (ITYPE.EQ.3)  THEN 
WRITE(NOUT, 1600) 

DO  50  I  «  1,  IPARAM 

READ(NINP,*)  ANOM(I) ,  AKLO(I) ,  AKHI(I) 

ALO(I)  =  ANOM(I)/AKLO(I) 

AHI(I)  =  ANQM(I)*AKHI(I) 

50  CONTINUE 

C 

C  ECHO  PARAMETER  INPUT  DATA  AND  COMPUTED  RANGE. 

C 

WRITE(NPRT, 1700) 

WRITE(NOUT, 1700) 

DO  60  I  *  1,  IPARAM 

WRITE(NPRT, 1800)  I,  ANOM(I) ,AKLO(I) ,AKHI(I) ,ALO(I) ,AHI(I) 
WRITE (NOUT, 1800)  I,  ANOM(I) , AKLO(I) ,AKHI(I) ,ALO(I) ,AHI(I) 
60  CONTINUE 

C 

C  ENTER  DATA  FOR  PARAMETERS  WHICH  ARE  UNIFORMLY  DISTRIBUTED. 

C 

ELSE  IF  (ITYPE.EQ.4)  THEN 
WRITE(NOUT, 1900) 

DO  70  I  *  i,  IPARAM 

READ(NINP ,*)  ANOM(I) ,  AKLO(I) ,  AKHI(I) 
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ALQ(I)  =  ANOM(I)/AKLO(I) 

AHI(I)  =  ANOM(I) *AKHI (I) 

70  CONTINUE 
C 

C  ECHO  PARAMETER  INPUT  DATA  AND  COMPUTED  RANGE. 

C 

HRITE(NPRT , 2000) 

WRITE(NOUT, 2000) 

DO  80  I  *  1,  IPARAM 

WRITE(NPRT,2100)  I,  ANOM(I) ,AKLO(I) ,AKHI(I) ,ALO(I) ,AHI(I) 
WRITE(N0UT,2100)  I,  ANOM(I) ,AKLO(I) ,AKHI(I) ,ALO(I) ,AHI(I) 

80  CONTINUE 
END  IF 
C 

C  CHECK  THE  VALUE  OF  'IOPT1  TO  DETERMINE  IF  THE  NUMBER  AND/OR 

C  VALUES  OF  THE  INDEPENDENT  PARAMETERS  USED  BY  THE  ALGORITHM 

C  NEED  TO  BE  ALTERED. 
r 

90  IF  ((I0PT.EQ.2) .OR. (I0PT.EQ.3) .OR. (IOPT. EQ. 4) .OR. (IOPT.Eq.8) .OR. 
1  (IOPT. Eq. 9). OR. (IOPT. EQ. 10). OR. (IOPT. Eq. 14))  GO  TO  110 
C 

C  SPECIFY  THE  NUMBER  'IXVAL'  OF  INDEPENDENT  PARAMETERS  WHICH 

C  APPEAR  IN  THE  ALGORITHM. 

C 

WRITE (NOUT, 2200) 

READ(NINP,*)  IXVAL 
C 

C  ENTER  THE  VALUES  OF  THE  'IXVAL'  INDEPENDENT  PARAMETERS. 

C 

WRITE (NOUT, 2300) 

READ(NINP,*)  (XVAL(I) ,  1*1,  IXVAL) 

C 

C  ECHO  THE  INDEPENDENT  PARAMETER  VALUES. 

C 

WRITE (NPRT, 2400) 

WRITE (NOUT, 2400) 

DO  100  I  =  1,  IXVAL 

WRITE (NPRT, 2500)  I,  XVAL(I) 

WRITE (NOUT, 2500)  I,  XVAL(I) 

100  CONTINUE 
110  CONTINUE 
RETURN 

1000  FORMAT(/ , '  ENTER  THE  NOMINAL  VALUES  OF  THE  PARAMETERS  ',/, 

1  '  WHICH  ARE  NORMALLY  DISTRIBUTED  AND  THE  ’,/, 

2  '  CORRESPONDING  K-FACTOR.' ,/) 

1100  FORMAT(/,’  PARAMETER  NOMINAL  K-FACTOR  95%  EXTREMES',/) 
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1200  F0RMAT(1X 
1300  FORMAT (/ , 
1 
2 

1400  FORMAT (/, 
1500  FORMAT (IX 
1600  FORMAT (/, 
1 
2 
3 

1700  FORMAT (/ , 
1 

1800  FORMAT (IX 

1  IX 
1900  FORMAT (/, 

1 

2 

2000  FORMAT (/, 
1 

2100  FORMAT ( IX 
1  IX 

2200  FORMAT (/ , 
1 

2300  FORMAT (/, 
2400  FORMAT (/, 
1 


13, 5X, 1PE10 .4, 2X.0PF10 .3,2X, 1PE10 . 4,2X, 1PE10 .4) 

ENTER  THE  NOMINAL  VALUES  OF  THE  PARAMETERS  ’  ,/, 

WHICH  ARE  LOGNORMALLY  DISTRIBUTED  AND  THE  ' ,/, 
CORRESPONDING  K-FACTOR.’,/) 

PARAMETER  NOMINAL  K-FACTOR  95*/.  EXTREMES ’, /) 

13 , 5X, 1PE10 . 4 , 2X , 0PF10 . 3 , 2X , 1PE10 . 4 , 2X , 1PE10 . 4) 

ENTER  THE  NOMINAL  VALUES  OF  THE  PARAMETERS  ’ ,/, 

WHICH  ARE  BETA  DISTRIBUTED  AND  THE  K-FACTORS , ' , / , 
(BOTH  K-LOW  AND  K-HIGH,  EVEN  IF  EQUAL,  MUST',/, 

BE  ENTERED) . ' ,/) 

PARAMETER  NOMINAL  K-FACTORS ' , 

EXTREMES',/) 

13 , 5X, 1PE10 . 4 , 2X, 0PF10 . 3 , 2X, 0PF10 . 3 , 2X, 1PE10 . 4 , 
1PE10.4) 

ENTER  THE  NOMINAL  VALUES  OF  THE  UNIFORMLY  ' , / , 
DISTRIBUTED  PARAMETERS  AND  THE  K-FACTOPS , ' , / , 

(ENTER  BOTH  K-LOW  AND  K-HIGH).',/) 

PARAMETER  NOMINAL  K-FACTORS', 

EXTREMES’,/) 

I3,5X, 1PE10 .4,2X,0PF10 .3,2X,0PF10 . 3,2X, 1PE10 .4, 
1PE10.4) 

ENTER  THE  NUMBER  OF  INDEPENDENT  VARIABLES  (IXVAL)',/, 
WHICH  APPEAR  IN  THE  ALGORITHM.’,/) 

ENTER  THE  VALUES  OF  THE  INDEPENDENT  VARIABLES.',/) 
INDEPENDENT  ’,/, 

PARAMETER  VALUE  ’,/) 


2500  F0RMAT(5X, 15 ,5X, 1PE12 .4) 
END 
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SUBROUTINE  ASSGNBETA(ANOM ,  ALO,  AHI,  AKLO,  AKHI ,  ASIG, 

1  Y,  AREA,  IBD) 

C 

C  SUBROUTINE  ASSGNBETA  PROMPTS  THE  USER  TO  SPECIFY  WHETHER 

C  TABULATED  OR  CALCULATED  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS 

C  WILL  BE  USED  TO  MODEL  PARAMETER  AND  ALGORITHMIC  UNCERTAINTIES. 

C  THE  CONTROL  PARAMETER  WHICH  DICTATES  WHICH  OF  THE  TWO  WILL  BE 
C  USED  IS  ’ ITABL ’  AND  IS  READ  BY  THE  SUBROUTINE.  THE  ACTION 
C  RESULTING  FROM  A  SPECIFIED  POSSIBLE  VALUE  OF  ITABL  IS  AS  FOLLOWS: 

C 

C  ITABL 

C 
C 
C 
C 

C  ITABL 

C 
C 
C 
C 

COMMON  /BLKOl/  IPARAM,  ITYPE,  ISIZE,  ISEED,  IALGO ,  IXVAL,  ITABL, 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

REAL*4  ANOM(*),  ALO(*) ,  AHI(*),  AKLO(*) ,  AKHI(*) 

REAL*4  ASIG(2,*) ,  Y(*) ,  AREA(*) 

INTEGERS  IBD(*) 

C 

C  INQUIRE  WHETHER  TABULATED  OR  CALCULATED  BETA  CUMULATIVE 

C  PROBABILITY  DISTRIBUTIONS  ARE  TO  BE  USED. 

C 

WRITE (NOUT, 1000) 

1000  FORMAT (/ , '  SELECT  WHETHER  TO  USE  TABULATED  OR  CALCULATED' ,/ , 

1  ’  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTION  VALUES.',/, 

2  '  1  -  USE  TABULATED  BETAS.',/, 

3  '  2  -  USE  CALCULATED  BETAS.’,/) 

READ (NINP,*)  ITABL 

C 

C  FOR  ITABL  *  1,  DIVERT  TO  SUBROUTINE  BETATABL  TO  ASSIGN  THE 

C  APPROPRIATE  TABULATED  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS. 

C 

IF  (ITABL. EQ.l)  THEN 

CALL  BETATABL (ANOM,  AKLO,  AKHI,  IBD) 

C 

C  FOR  ITABL  *  2,  DIVERT  TO  SUBROUTINE  BETACALC  TO  CALCULATE 

C  THE  REQUIRED  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS. 

C 


=  1:  SUBROUTINE  BETATABL  IS  CALLED  IN  ORDER  TO  ASSIGN 

THE  APPROPRIATE  TABULATED  BETA  CUMULATIVE  PROBABILITY 
DISTRIBUTION  TO  EACH  PARAMETER  AND  ALGORITHM  WITH 
UNCERTAINTY . 

=  2:  SUBROUTINE  BETACALC  IS  CALLED  TO  CALCULATE  THE  SPECIFIC 
BETA  CUMULATIVE  PROBABILITY  DISTRIBUTION  WHICH  SHOULD 
BE  ASSOCIATED  WITH  EACH  PARAMETER  AND  ALGORITHM  WITH 
UNCERTAINTY . 
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ELSE  IF  (ITABL.EQ .2)  THEN 

CALL  BETACALC(ANOM,  ALQ,  AHI ,  ASIG,  Y,  AREA,  IBD) 
END  IF 
RETURN 
END 
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SUBROUTINE  BETATABL (ANOM ,  AKLO ,  AKHI ,  IBD) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  BETATABL  USES  THE  USER  SPECIFIED  PARAMETER  HIGH-SIDE 
AND  LOW-SIDE  K-FACTOR  UNCERTAINTIES  TO  DETERMINE  THE  TABULATED 
BETA  CUMULATIVE  PROBABILITY  DISTRIBUTION  WHICH  WILL  BE  USED  TO 
MODEL  THE  UNCERTAINTY  ASSOCIATED  WITH  EACH  PARAMETER.  THE 
SPECIFIED  K-FACTORS  FOR  A  PARAMETER  ARE  FIRST  USED  TO  CALCULATE 
THE  'MEAN'  (TEMP)  OF  THE  ASSOCIATED  DISTRIBUTION: 

TEMP  =  (1.  -  1./ AKLO) /(AKHI  -  l./AKLO) 

THIS  'MEAN'  IS  THEN  COMPARED  TO  THE  'MEAN'  (ll-TH  ELEMENT)  OF  EACH 
TABULATED  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTION.  THE  LATTER 
WHICH  HAS  A  'MEAN'  CLOSEST  TO  THAT  DETERMINED  WITH  THE  SPECIFIED 
PARAMETER  K-FACTORS  IS  THEN  USED  TO  MODEL  THE  UNCERTAINTY 
ASSOCIATED  WITH  THAT  PARAMETER. 

THE  MAIN  ARRAYS  USED  BY  THE  SUBROUTINE  ARE: 

BD  -  ARRAY  CONTAINING  81  TABULATED  BETA  CUMULATIVE  PROBABILITY 
DISTRIBUTIONS.  EACH  DISTRIBUTION  CONSISTS  OF  21  VALUES 
WHICH  CORRESPOND  TO  THE  FRACTION  OF  THE  UNIT  INTERVAL 
(FROM  0  TO  1)  AT  WHICH  5’/.  INCREMENTS  IN  THE  CUMULATIVE 
PROBABILITY  ARE  REACHED.  FOR  INSTANCE,  IF  THE  FIRST  THREE 
ENTRIES  IN  THE  SET  OF  DATA  FOR  A  DISTRIBUTION  ARE  O.O,  0.113, 
AND  0.187,  THEN  THE  CUMULATIVE  PROBABILITY  ASSOCIATED  WITH 
THE  ENTRIES  IS  0.0,  0.05,  AND  0.10,  RESPECTIVELY. 

IBD  -  INTEGER  ARRAY  WHOSE  I-TH  ENTRY  IS  THE  NUMBER  OF  THE  TABULATED 
BETA  CUMULATIVE  PROBABILITY  DISTRIBUTION  WHICH  IS  IDENTIFIED 
AS  HAVING  A  MEDIAN  NEAREST  IN  VALUE  TO  THE  COMPUTED  MEDIAN 
OF  THE  I-TH  PARAMETER  WHICH  HAS  BETA  DISTRIBUTED  UNCERTAINTY. 
THE  TABULATED  DISTRIBUTION  THUS  IDENTIFIED  IS  USED  TO  MODEL 
THE  UNCERTAINTY  ASSOCIATED  WITH  THE  PARAMETER. 

ANOM  -  ARRAY  OF  PARAMETER  NOMINAL  VALUES. 

AKLO  -  ARRAY  OF  PARAMETER  LOW-SIDE  K-FACTOR  UNCERTAINTIES. 

AKHI  -  ARRAY  OF  PARAMETER  HIGH-SIDE  K-FACTOR  UNCERTAINTIES. 

COMMON  /BLK01/  IPARAM,  I TYPE,  ISIZE,  ISEED,  IALGO ,  IXVAL,  ITABL , 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

COMMON  /BLK02/  BD(21,81) 

REAL*4  ANOM(*) ,  AKLQ(*) ,  AKHI (*) 

REAL*4  DBA(21 ,41) ,  DBB(21,40) 

INTEGERS  IBD(*) 
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DATA  DBA/O . , . 015 , . 029 , . 044, 


.  059 , . 075 , . 091 , . 109 , . 127 , . 146 , . 167 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


, .189, .213, .239,. 269,. 302, .340, .386, .446, .535,1.  , 
0. , .016, .031,-047, .062, .079, .096, . 113, . 132, . 151, . 172 
, .195, .219, .245,. 275,. 308,. 346, .392 , .452 , . 540 , 1 . , 
0.  ,  .022,-038, .054, .070 , . 087 , . 104, . 122, . 141, . 161 , . 182 
, .205, .229, .255, .285, .318, .356, .402, .461 , . 548, 1 . , 
0. , .020, .037, .053, .070, .088, . 105, . 124, . 143, . 163, . 185 
, .208, .232, .259,. 289,. 322, .360,-406, .465, .552,1 . , 
0. , .024, .042, .060, .077, . 095 , . 113, . 132, . 151, . 172 , . 194 
, .217, .241, .268, .298, .331, .369, .415, .474, . 560 , 1 . , 
0. ,.024, .044,. 062, .081, .099, .118, .137, .157, .178,. 200 
, .224, .249, .276,. 305, .339, .377, .423, .481, .566,1. , 
0. , .029, .050, .069, .088, . 107, . 126, . 146, . 166, . 188,  .210 
, .233, .259, .286, .316, .349, .387, .432, .490, .575, 1 .  , 
0. , .030, .052, .072, .092, .112, .132, .152, .173, .194,. 217 
, .241, .266, .293,-323, .357, .395 , .440 , .498 , . 581 , 1 .  , 
0. , .035, .058, .080, .100, .120, .141, .162, .183, .205,-227 
, .251, .277, .304,. 334, .368, .406, .450,-507, .590,1. , 
0. , .039, .064, .087, . 108, . 129, . 150, . 171 ,. 193, . 215, . 238 
, .263, .288, .316, .346, .379, .417, .461 , . 518, . 599 , 1 .  , 
0. , .043, .070, .094,. 117, .138, .160,. 182, .204, .227,. 250 
, .274, .300, .328,. 358,. 391, .429,-473, .529, .609,1.  , 
0. , .046, .074, .098, .121, .143, .165, .187, .210, .232, .256 
, .281, .306, .334,. 364,. 397, .435, .479, .534, .614,1 .  , 
0. , .050, .078, .103, .127, .150, .172, .194, .217, .240, .263 
, .288, .314, .342, .372, .405, .442, .486 , .541 , . 620 , 1 . , 
0. , .053, .083, .109, .133, .156, .178, .201, .224, .247,. 271 
, .295, .321, .349, .379,-412, .449,-493, .547, .625, 1 .  , 
0. , .057, .087,. 114, .138, .162, .185,. 208, .230, .254, .278 
, .303, .329, .356, .386, .419, .456, .499, .553, . 631 , 1 . , 

0 . , .061 , . 093 , . 120 , . 145 , . 169 , . 192 , . 215, . 239 , . 262 , . 286 
, .311, .337, .365, .395, .427, .464, .507, .560 , . 637 , 1 .  , 
0. , .065, .098, .126, .152, .176, .200, .223,-246, .270, .294 
, .319, .345, .373,-403, .435, .472, .514, .567, .643,1.  , 

0 . , . 070 , . 105 , . 134 , . 160 , . 184 , . 208 , . 232 , . 255 , . 279 , . 303 
, .329, .355, .382, .412, .444, .481 , .523, . 575 , . 650 , 1 .  , 
0. , .076, .111, .141, .168, . 193, . 217 , . 241 , . 264, . 288 , . 313 
, .338, .364, .391, .421, .453, .489 , .531 , . 583 , . 656 , 1 .  , 
0. , .082, .119,-149, .176, .202, .226, .250, .274, .298,. 323 
, .348, .374, .401, .431, .463, .498, .540 , . 591 , . 664 , 1 .  , 
0. , .089, . 127, .158, .186, . 212, . 237, . 261 , . 285, .309, .334 
, .359, .385, .412, .441, .473, .508, .549 , . 600 , . 671 , 1 .  , 
0. , .092, .131, .162, .190, .216, .241, .266, .290, .314,-338 
, .363, .389, .417,. 446,. 478, .513, . 554 , . 604, . 675 , 1 .  , 

0 . ,  . 096 , . 136 , . 168 , . 196 , . 223 , . 248 , . 273, . 297, . 321 , . 345 
, .370, .396,.424, .453, .484, .519 , .560 , . 609 , . 680 , 1 .  , 
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.263, .324, .367, .403, .435, 
.589, .613, .638, .663, .690, 
.269, .330, .374, .411,-443, 
.598, .622, .647,. 672,. 699, 
.274,. 337,. 381,. 418,. 451, 
.607, .632, .656, .681, .708, 
.279, .343, .388, .426, .458, 
.616, .640, .665, .690, .716, 
.285, .349, .395, .433, .466, 
.625, .649, .673, .699, .725, 
.290, .355, .402, .440, .473, 
.633, .657, .681,. 706,. 733, 
.294, .360, .407, .446,-480, 
.640, .664, .689,. 714,. 740, 
.299, .366, .413, .452, .486, 
.647, .671, .696, .721, .747, 
.304, .371, .420, .459, .493, 
.654, .679, .703, .728, .754, 
.308, .376, .425, .465, .499, 
.661, .685, .710, .734, .760, 
.313, .382, .431, .471, .505, 
.668, .692, .716, .741, .766, 
.317, .386, .435, .476, .511, 
.673, .698, .722,. 746,. 772, 
.320,. 391,. 440,. 481,. 516, 
.679, .703, .727,. 752,. 777, 
.325,-396, .446, .487, .522, 
.686, .710, .734, .759, .784, 
.329, .400, .451, .492, .527, 
.691, .715, .739, .763, .788, 
.336, .409, .460, .502, .537, 
.702, .726, .750, .774, .798, 
.344, .417, .469, .511, .547, 
.712, .736, .759, .783, .807, 
.350, .425, .477, .519,-556, 
.721, .745, .768, .792, .816, 
.357, .433, .486, .528, .565, 
.730, .754, .777,. 800,. 824, 
.363, .440, .493, .536, .573, 
.738, .761, .785, .808, .831, 
.369,. 447,. 501,. 544,. 581, 
.746, .770, .792, .815, .838, 
.375,. 453,. 507,. 551,. 588, 
.753, .776, .799,. 822,. 844, 
.380, .459, .514, .558, .595, 
.760, .783, .806, .828, .850, 


.464, .491,-516, .541, .565 
.718, .750, .787, .836,1. , 
.472, .499, .525, .550,. 574 
.727, .758, .795, .843,1. , 
.480, .508, .534, .559,. 583 
.736, .767, .803, .851,1. , 
.488, .516, .542, .567,. 592 
.744, .775, .811, .857,1. , 
.496, .524, .550, .576..600 
.753, .783, .819, .864,1. , 
.503, .531, .558, .583,. 608 
.760, .791, .826, .871,1. , 
.510, .538, .565, .590, .615 
.767,. 797,. 832,. 876,1., 
.517, .545, .572, .598,-623 
.774, .804, .838, .881,1. , 
.524, .552,.579, .605, .630 
.781, .810, .844, .887,1. , 
.530, .559, .586, .612, .637 
.787, .816, .849, .891,1. , 
.536, .565, .592, .618, .643 
.793, .822, .855, .896,1. , 
.542, .571, .598, .624,. 649 
.798, .827, .859, .900,1. , 
.547, .576, .604, .630,. 655 
.804, .832, .864, .904,1. , 
.554, .583, .611, .637,. 662 
.810, .838, .869, .908,1. , 
.559, .588, .615, .641, .666 
.814, .842, .873, .911,1. , 
.569, .599, .626, .652, .677 
.824, .851, .881, .918,1. , 
.579, .609, .636, .662, .687 
.832, .859, .889, .924,1. , 
.588, .618, .645, .671, .697 
.840, .866, .895, .930,1. , 
.597, .627, .655, .681, .706 
.848, .874, .902, .935,1. , 
.605,-635, .663, .689, .714 
.855, .880, .907, .939,1 . , 
.614, .644, .671, .697, .722 
.862,  .886, .913, .943,1. , 
.621, .651, .679, .705, .729 
.867, .891, .917, .947,1. , 
.628, .658, .686, .712,-737 
.873, .897, .922, .950,1. , 
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*  0. , .386, .466, .521, .565, .603, .636, .666, .694, .719,. 744 

*  , .768, .790, .813, .835, .857, .879, .902,-926, .954,1.  , 

*  0. , .391, .471, .527, .571, .609, .642, .672, .700, .726,. 750 

*  , .773, .796, .818, .840, .862, .883, .906, . 930 , .957 , 1 .  , 

*  0. , .401, .483, .539, .584, .622, .655, .685, .712, .738, .762 

*  , .785, .808, .829, .851,. 872, .893, .914, .937, .962,1.  , 

*  0. , .410, .493, .550, .595, .633, .666, .696, .724, .749,-773 

*  , .796, .818, .839, .860, .881, .901, .922, .943, .966,1.  , 

*  0. , .419, .503, .560, .605, .643, .677, .707, .734, .759,-783 

*  , .806, .827, .848, .868, .888, .908, .928, .948, .970,1.  , 

*  0. , .426, .511, .568, .614, .652, .685, .715, .743, .768, .791 

*  , .814, .835, .855,. 875,. 895, .914, .933, .952, .973, 1 .  , 

*  0. , .434, .519, .577, . 623, . 661 , . 695 , .724, .751, .776,-800 

*  , .822, .843, .863, .882,. 901, .920, .938, . 956 , . 976 , 1 .  , 

*  0. , .441,-528, .586, .632, .671, .704, .734, .761, .785, .808 

*  , .830, .851, .870, .889, .908, .925, .943, .960, .979,1.  , 

*  0. , .448, .535,-594, .640, .678, .712, .741, .768, .793, .815 

*  , .837, .857, .876, .895, .913, .930, .947, .964, .981,1.  , 

*  0. , .453, .541, .600, .646, .684, .718, .747, .774, .798,. 821 

*  , .842, .862,. 881,. 899,. 917, .933, .950, . 966 , . 982 , 1 .  , 

*  0. , .460, .548, .608, .654, .692, .726, .755, .781, .806, .828 

*  , .849, .869, .887, .90S,. 922, .938, .954, .969, .984,1. , 

*  0. , .465, .554, .614, .660, .699, .732, .761, .787, .811, .833 

*  , .854, .873, .892,. 909,. 925,. 941, .956, .971, .985,1./ 
EQUIVALENCE  (BD(1 ,1) ,DBA(1 , 1)) 

EQUIVALENCE  (BD(1 ,42) ,DBB(1 , 1) ) 

BDMAX=20 . 0 
C 

C  LOOP  OVER  THE  PARAMETERS  WITH  UNCERTAINTY  TO  BE  MODELED  WITH 

C  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS. 

C 

DO  20  I  »  1,  IPARAM 
C 

C  EVALUATE  THE  'MEAN'  =  TEMP  OF  THE  PARAMETER  AS  DETERMINED  FROM 

C  ITS  SPECIFIED  K-FACTOR  UNCERTAINTIES. 

C 

TEMP  =  (1.0-1 .0/AKL0(I))/ (AKHI(I)-l . 0/AKL0(I)) 

C 

C  LOOPING  OVER  ALL  THE  TABULATED  BETA  CUMULATIVE  DISTRIBUTIONS, 

C  COMPARE  THE  PARAMETER  COMPUTED  'MEAN'  TO  THOSE  (TEMPI  AND  TEMP2) 

C  OF  THE  TABULATED  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTIONS.  THE 

C  'MEAN’  OF  THE  J-TH  TABULATED  DISTRIBUTION  IS  THE  11-TH  ELEMENT  IN 

C  THE  ARRAY  BD(11,J).  IF  THE  J-TH  TABULATED  DISTRIBUTION  HAS  A 

C  'MEAN'  CLOSEST  TO  THE  COMPUTED  'MEAN'  OF  THE  I-TH  PARAMETER,  ASSIGN 

C  THAT  DISTRIBUTION  TO  THE  I-TH  PARAMETER  (IBD(I)  =  J) . 

C 
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TEMPI  =  0 .5*(BD(11 , l)+BD(ll ,2) ) 

TEMP2  =  0 . 5*(BD(ll ,80)+BD(ll,81) ) 

IF  (TEMP. LE. TEMPI)  THEN 
IBD(I)  *  1 

ELSE  IF  (TEMP . GT . TEMP2 )  THEN 
IBD(I)  *  81 
ELSE 

DO  10  J  *  2,  80 

TEMPI  =  0.5*(BD(ll,J)+BD(ll,J-l)) 

TEMP2  *  0 .5*(BD(11, J)+BD(ll , J+l) ) 

IF  ( (TEMP. GT. TEMPI) .AND. (TEMP. LE.TEMP2))  IBD(I)  =  J 
10  CONTINUE 

END  IF 
20  CONTINUE 
C 

C  OUTPUT  THE  NUMBER  OF  THE  CUMULATIVE  BETA  PROBABILITY 

C  DISTRIBUTION  ASSIGNED  TO  EACH  PARAMETER  WITH  UNCERTAINTY. 

C 

WRITE (NPRT, 1000) 

WRITE  (NOUT, 1000) 

DO  30  I  *  1,  IPARAM 

WRITE(NPRT.llOO)  I,  ANOM(I) ,  AKLO(I),  AKHI(I) ,  IBD(I) 

WRITE (NOUT, 1100)  I,  ANOM(I) ,  AKLO(I),  AKHI(I),  IBD(I) 

30  CONTINUE 
RETURN 

1000  FORMAT(/ , ’  SELECTED  TABULATED  BETA  CUMULATIVE  PDF  » ,/, 

1  »  PARAMETER  NOMINAL  UNCERTAINTY  FACTORS  BETA  ' , / , 

2  »  NUMBER  VALUE  K-LOW  K-HIGH  NUMBER’,/) 

1100  F0RMAT(3X, 15 , 2X , 1PE10 . 4 , 2X , 0PF10 . 3 , 2X, 0PF10 . 3 , 5X , 15) 

END 
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O  V 


SUBROUTINE  BETACALC (ANOM ,  ALO,  AHI ,  ASIG,  Y,  AREA,  IBD) 

C 

C  SUBROUTINE  BETACALC  IS  USED  TO  COMPUTE  THE  BETA  CUMULATIVE 

C  PROBABILITY  DISTRIBUTION  TO  BE  USED  TO  MODEL  THE  UNCERTAINTY 

C  ASSOCIATED  WITH  EACH  PARAMETER  WITH  BETA  DISTRIBUTED  UNCERTAINTY. 

C  FIRST  THE  NOMINAL,  LOW,  AND  HIGH  VALUES  OF  THE  PARAMETER  ARE 

C  USED  TO  COMPUTE  THE  'MEAN'  M  OF  THE  DESIRED  DISTRIBUTION: 

C 

C  M  =  (ANOM  -  ALO)/ (AHI  -  ALO) 

C 

C  THIS  'MEAN'  IS  THEN  USED  TO  COMPUTE  THE  POSSIBLE  MAXIMUM 
STANDARD  DEVIATIONS  (SIG1  AND  SIG2)  FOR  THE  PARAMETER: 

C  SIG1  =  M*SqRT((l  -  M)/(l  +  M)) 

C  SIG2  «  (1  -  M)*SQRT(M/ (2  -  M)) 

C 

C  THE  USER  IS  THEN  REQUESTED  TO  SPECIFY  THE  STANDARD  DEVIATION  SIG 

C  FOR  THE  PARAMETER  WHICH  MUST  SATISFY  THE  FOLLOWING: 

C 

C  SIG  LESS  THAN  OR  EQUAL  TO  MIN(SIG1,  SIG2) 

C 

C  THE  SPECIFIED  STANDARD  DEVIATION  FOR  THE  PARAMETER  IS  NEXT  USED 

C  TO  COMPUTE  A  ’NORMALIZED'  STANDARD  DEVIATION  SIGMA: 

C 

C  SIGMA  =  SIG/ (AHI  -  ALO) 

C 

C  AND  THE  TWO  PARAMETERS  OF  THE  BETA  DISTRIBUTION  (ALPHA  AND  BETA) 

C  ARE  THEN  COMPUTED: 

C 

C  ALPHA  =  M*[{(M  -  M*M)/SIGMA**2}  -  1] 

C  BETA  =  (1  -  M)*[{(M  -  M*M)/SIGMA**2>  -  1] 

C 

C  FUNCTION  ROUTINE  FX  IS  NEXT  USED  TO  EVALUATE  THE  BETA  DISTRIBUTION 

C  CHARACTERIZED  BY  THE  PARAMETERS  ALPHA  AND  BETA,  AND  THE  ASSOCIATED 

C  CUMULATIVE  PROBABILITY  DISTRIBUTION  (AT  5’/.  AREA  INCREMENTS)  IS 

C  DETERMINED  AND  SAVED  FOR  USE  IN  MODELING  THE  PARAMETER  UNCERTAINTY. 

C 

COMMON  /BLKOl/  IPARAM,  I TYPE,  ISIZE,  I SEED,  I ALGO ,  IXVAL,  ITABL, 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

COMMON  /BLK02/  BD(21,81) 

REAL*4  ANOM(*) ,  ALO(*),  AHI(*) 

REAL*4  ASIG(2 , *) ,  Y(*)  ,  AREA(*) 

INTEGER*4  IBD(*) 

C 

C  DEFINE  THE  INCREMENT  USED  IN  EVALUATING  THE  BETA  DISTRIBUTION  AND 
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C  CALCULATING  THE  AREA  INCREMENTS  UNDER  THE  DISTRIBUTION. 

C 

DELINC  =  l.O/NINC 
C 

C  COMPUTE  THE  MAXIMUM  ALLOWABLE  STANDARD  DEVIATION  FOR 
C  THE  PARAMETER.  LOOP  OVER  ALL  PARAMETERS  WITH  UNCERTAINTY. 

C 

DO  20  I  -  1.  IPARAM 

TEM  *  (ANOM(I)  -  ALO(I))/(AHI(I)  -  ALO(I)) 

ASIG(l.I)  *  TEM*SQRT((1.0  -  TEM)/(1.0  +  TEM)) 

ASIG(1 , I)  =  (AHI(I)  -  ALO (I) ) *ASIG (1,1) 

ASIG(2 , I)  *  (1.0  -  TEM)*SQRT(TEM/ (2 . 0  -  TEM)) 

ASIG(2,I)  -  (AHI(I)  -  AL0(I))*ASIG(2,I) 

20  CONTINUE 

DO  30  I  *  1,  IPARAM 

WRITE(NPRT, 1000)  I,  ASIG(1,I),  ASIG(2,I) 

WRITE (NOUT, 1000)  I,  ASIG(1,I),  ASIG(2,I) 

30  CONTINUE 
C 

C  INTERACTIVELY  READ  THE  PARAMETER  SPECIFIED  STANDARD  DEVIATION. 

C 

WRITE (NPRT, 1100) 

WRITE(NQUT, 1100) 

READ(NINP,*)  (ASIG(l , I) ,  1*1,  IPARAM) 

C 

C  ECHO  THE  INPUT  VALUES. 

C 

WRITE (NPRT, 1200) 

WRITE (NOUT, 1200) 

DO  40  I  *  1,  IPARAM 

WRITE (NPRT, 1300)  I,  ANOM(I) ,  ALO(I),  AHI(I),  ASIG(1,I) 

WRITE (NOUT, 1300)  I,  ANOM(I) ,  ALO(I),  AHI(I),  ASIG(1,I) 

40  CONTINUE 

WRITE (NPRT, 1400) 

WRITE (NOUT, 1400) 

C 

C  LOOP  OVER  PARAMETERS  WITH  UNCERTAINTY. 

C 

DO  90  I  *  1,  IPARAM 
C 

C  COMPUTE  THE  NORMALIZED  PARAMETER  STANDARD  DEVIATION,  PARAMETER 
C  NORMALIZED  MEAN,  AND  THE  ALPHA  AND  BETA  PARAMETERS  FOR  THE 

C  APPROPRIATE  BETA  DISTRIBUTION. 

C 

ASIG(1 , I)  =  ASIG(l , I) / (AHI (I)  -  ALO(I) ) 

TEM  *  (ANOM(I)  -  ALO (I))/ (AHI (I)  -  ALO (I)) 
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TEM1  *  (TEM  -  TEM*TEM)/ (ASIG(1 , I)*ASIG(1 , I) )  -  1.0 

ALPHA  =  TEM*TEM1 

BETA  *  (1.0  -  TEM) *TEM1 

OUTPUT  THE  COMPUTED  ALPHA  AND  BETA  PARAMETERS. 

WRITE (NPRT, 1500)  I,  ALPHA,  BETA 
WRITE (NOUT, 1500)  I,  ALPHA,  BETA 

COMPUTE  THE  BETA  DISTRIBUTION  WITH  PARAMETER  VALUES 
ALPHA  AND  BETA.  USE  THE  1AAPEZ0IDAL  RULE  TO  COMPUTE 
THE  TOTAL  AREA  (SUMNORM)  UNDER  THE  COMPUTED  BETA 
DISTRIBUTION. 


Y(l)  =  0.0 
SUMNORM  =0.0 
DO  50  J  =  1,  NINC 

Y(J+l)  =  Y(J)  +  DELINC 
FX1  =  FX (ALPHA,  BETA,  Y(J)) 

FX2  =  FX  (ALPHA,  BETA,  Y(J+D) 

SUMNORM  =  SUMNORM  +  0.5*(FX1  +  FX2)*(Y(J+1)  -  Y(J)) 
CONTINUE 

COMPUTE  THE  NORMALIZED  BETA  DISTRIBUTION. 

DO  60  J  =  1,  NINC+1 

FX1  =  FX (ALPHA,  BETA,  Y(J)) 

Y(J)  =  FX1/SUMN0RM 
CONTINUE 

EVALUATE  AREA  INCREMENTS  FOR  THE  AREA  UNDER  THE  NORMALIZED 
BETA  DISTRIBUTION. 

AREA(l)  =  0 . 5*DELINC* (Y ( 1 )  +  Y(2)) 

DO  70  J  =  2,  NINC 

AREA(J)  =  AREA(J-l)  +  0.5*DELINC*(Y(J)  +  Y(J+l)) 
CONTINUE 

DETERMINE  THE  ABSCISSA  VALUES  FOR  THE  57.  AREA  INCREMENTS 
(HENCE  THE  CUMULATIVE  PROBABILITY  DISTRIBUTION)  FOR  THE 
COMPUTED  NORMALIZED  BETA  DISTRIBUTION.  SAVE  THE  COMPUTED 
VALUES  IN  THE  ARRAY  BD(J,I)  (J  =  1,  21)  FOR  THE  I-TH 
PARAMETER. 

BD(1 , I)  =  0.0 
BD(21 , I)  =  1.0 


IBD(I)  *  I 
DO  90  J  *  2,  20 

AAREA  *  0.05*(J  -  1) 

DO  80  K  *  1,  NINC-1 
IF  (AAREA. GT. AREA (K+l))  THEN 
GO  TO  80 

ELSE  IF  ((AAREA. GT.AREA(K)) .AND. (AAREA. LE. AREA (K+l)))  THEN 
FRAC  *  (AAREA  -  AREA(K) ) / (AREA(K+l)  -  AREA(K)) 

BD(J,I)  =  K*DELINC  +  FRAC*DELINC 
GO  TO  90 

ELSE  IF  ( AAREA. LE. AREA (K))  THEN 
GO  TO  90 
ENDIF 

80  CONTINUE 

90  CONTINUE 
C 

C  OUTPUT  THE  BETA  CUMULATIVE  PROBABILITY  DISTRIBUTION  COMPUTED 
C  FOR  EACH  PARAMETER  WITH  UNCERTAINTY. 

C 

DO  110  I  «  1,  IP ARAM 
WRITE(NPRT, 1600)  I 
WRITE (NOUT, 1600)  I 
DO  100  J  *  1,  10 

FRAC1  *  (J  -  1)*0.05 
FRAC2  »  (J  +  10) *0.05 

WRITE (NPRT, 1700)  FRAC1,  BD(J,I),  FRAC2 ,  BD(J+11,I) 

WRITE (NOUT, 1700)  FRAC1 ,  BD(J,I),  FRAC2 ,  BD(J+11,I) 

100  CONTINUE 

WRITE (NPRT ,1700)  0.50,  BD(11,I) 

WRITE (NOUT, 1700)  0.50,  BD(11,I) 

110  CONTINUE 
RETURN 

1000  FORMAT ( ’  FOR  PARAMETER’ ,13, '  THE  STAN.  DEV.  MUST  BE  <’, 

1  '  MIN(’ ,1PE10.4, ' , ’ , 1PE10.4, '). ’) 

1100  FORMAT(/, '  ENTER  THE  STANDARD  DEVIATIONS  OF  THE  BETA',/, 

1  ’  DISTRIBUTED  PARAMETERS.’,/) 

1200  FORMAT(/,’  PARAMETER  NOMINAL  MINIMUM  MAXIMUM 
1  ’  STAN.  DEV.’,/) 

1300  FuRMAT(I5, 5X , 1PE10 . 4,2X , 1PE10 . 4, 2X , 1PE10 . 4, 2X , 1PE10 .4) 

1400  FORMAT(/,’  PARAMETER  ALPHA  BETA',/) 

1500  F0RMAT(3X, 15 ,5X,F7 .3, 5X ,F7. 3) 

1600  FORMAT (/ , ’  FOR  PARAMETER’ ,13, ’  THE  CUMULATIVE  BETA’,/, 

1  ’  DISTRIBUTION  VALUES  AT  5%  AREA  INCREMETJTS ' ,  /  , 

2  ’  ARE  AS  FOLLOWS: ’  ,/,/, 

3  '  AREA  AREA’,/, 

4  ’  FRACTION  VALUE  FRACTION  VALUE’,/) 
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1700  FORMAT ( 5X , F5 . 2 , 5X , F7 . 3 , 5X ,F5 . 2 , 5X , F7 . 3) 
END 
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SUBROUTINE  SOLVE(ANOM,  ALO,  AHI,  AKLO,  AKHI,  XVAL,  A,  IBD,  RESULT) 

SUBROUTINE  SOLVE  EVALUATES  THE  USER  ALGORITHM  ISIZE  TIMES  IN 
ORDER  TO  GENERATE  A  ALGORITHM  RESULTS  DISTRIBUTION  FROM  WHICH 
THE  OVERALL  COMBINED  UNCERTAINTY  TO  BE  ASSOCIATED  WITH  THE  OUTPUT 
OF  THE  ALGORITHM  IS  DETERMINED.  THE  ALGORITHM  IS  EVALUATED  IN  A 
TWO  STEP  PROCESS: 

1.  THE  PARAMETER  VALUES  TO  BE  USED  IN  EVALUATION  OF  THE  ALGORITHM 
ARE  RANDOMIZED.  A  NEW  RANDOM  NUMBER  (RN)  IS  USED  EACH  TIME 

TO  GENERATE  A  NEW  RANDOMIZED  VALUE  FOR  EACH  PARAMETER  WITH 
UNCERTAINTY  IN  THE  ALGORITHM.  EACH  RANDOM  NUMBER  USED  FOR 
OBTAINING  VALUES  FOR  PARAMETERS  WITH  NORMALLY  OR  LOGNORMALLY 
DISTRIBUTED  UNCERTAINTY  IS  NORMALLY  DISTRIBUTED  AND  OBTAINED 
BY  CALLING  SUBROUTINE  FNRN.  THE  RANDOM  NUMBERS  USED  FOR 
PARAMETERS  WHICH  ARE  BETA  OR  UNIFORMLY  DISTRIBUTED  ARE  OBTAINED 
BY  CALLING  SUBROUTINE  UPR1  AND  ARE  UNIFORMLY  DISTRIBUTED.  A  RANDO 
VALUE  FOR  A  PARAMETER  IS  A  FUNCTION  OF  THE  PARAMETER  NOMINAL 
VALUE,  THE  UNCERTAINTY  SPECIFICATION,  AND  THE  CURRENT  RANDOM 
NUMBER  VALUE.  IF  NECESSARY,  NOMINAL  VALUES  OF  PARAMETERS 
WHICH  ARE  GENERATED  FROM  SIMPLE  ALGORITHMS  ARE  FIRST  OBTAINED 
BY  CALLING  SUBROUTINE  AUXALGO.  THIS  WILL  ONLY  BE  DONE  IF  IXALG  = 

2.  THE  ALGORITHM,  WHICH  IS  CODED  IN  THE  USER  SUPPLIED  SUBROUTINE  TF, 
IS  THEN  EVALUATED  USING  THE  CURRENT  SET  OF  RANDOMIZED  PARAMETER 
VALUES.  THE  RESULT  OF  THE  EVALUATION  IS  RETURNED  TO  SUBROUTINE 
SOLVE  AND,  IF  NECESSARY  (IALGO  =  1),  ALGORITHMIC  UNCERTAINTY  IS 
APPLIED  TO  THE  RESULT.  THE  FINAL  COMPUTED  VALUE  FOR  THE  ALGORITHM 
IS  THEN  STORED  IN  ARRAY  RESULT. 

THE  ABOVE  TWO  STEP  PROCESS  IS  REPEATED  THE  USER  SPECIFIED  ISIZE  TIMES. 
THE  MAIN  VARIABLES  AND  ARRAYS  USED  BY  THE  SUBROUTINE  ARE  AS  FOLLOWS: 

ANOM  -  ARRAY  OF  NOMINAL  VALUES  OF  PARAMETERS  WITH  UNCERTAINTY. 

AKLO  -  ARRAY  OF  PARAMETER  LOW-SIDE  K-FACTOR  UNCERTAINTIES. 

AKHI  -  ARRAY  OF  PARAMETER  HIGH-SIDE  K-FACTOR  UNCERTAINTIES. 

ALO  -  ARRAY  OF  PARAMETER  COMPUTED  LOW-SIDE  EXTREME  VALUES. 

AHI  -  ARRAY  OF  PARAMETER  COMPUTED  HIGH-SIDE  EXTREME  VALUES. 

A  -  TEMPORARY  ARRAY  OF  RANDOMIZED  PARAMETER  VALUES. 

XVAL  -  ARRAY  OF  INDEPENDENT  PARAMETER  VALUES  (NO  UNCERTAINTY) . 

BD  -  ARRAY  OF  TABULATED  OR  CALCULATED  BETA  CUMULATIVE 

PROBABILITY  DISTRIBUTIONS . 

IBD  -  ARRAY  WHOSE  ENTRIES  IDENTIFY  THE  DISTRIBUTION  NUMBER 
ASSOCIATED  WITH  EACH  PARAMETER  WITH  BETA  DISTRIBUTED 
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UNCERTAINTY. 

RESULT  -  ARRAY  OF  EVALUATED  ALGORITHM  RESULTS. 

RN  -  VARIABLE  WHOSE  VALUE  IS  THE  CURRENT  NORMALLY  OR 
UNIFORMLY  DISTRIBUTED  RANDOM  NUMBER. 

COMMON  /BLKOl/  IPARAM,  ITYPE,  ISIZE,  ISEED,  IALGO ,  IXVAL,  ITABL, 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

COMMON  /BLK02/  BD(21,81) 

REAL*4  ANOM(*) ,  ALO(*),  AHI(*),  AKLO(*) ,  AKHI(*) 

REAL*4  XVAL(*) ,  A(*),  RESULT (*) 

INTEGER*4  IBD(*) 

BDMAX  =20.0 

INITIALIZE  THE  SEED  INTEGERS  FOR  THE  RANDOM  NUMBER  GENERATORS. 

NUMUR  =  ISEED 
NUMNR  =  ISEED 

SET  MAXIMUM  LOOP  COUNTER  VALUE  IPARA  DEPENDING  ON  WHETHER  OR  NOT 
THE  ALGORITHM  TO  BE  EVALUATED  HAS  OVERALL  ALGORITHMIC  UNCERTAINTY. 

IF  (IALGO .Eq. 1)  THEN 
IPARA  =  IPARAM  -  1 
ELSE 

IPARA  =  IPARAM 
END  IF 

OBTAIN  PARAMETER  NOMINAL  VALUES  FROM  SUBROUTINE  AUXALGO  IF  REQUIRED 
AND  THEN  COMPUTE  ASSOCIATED  LOW-SIDE  AND  HIGH-SIDE  VALUES. 

IF  ( IXALG. EQ.l)  CALL  AUXALGO (ANOM,  XVAL) 

IF  (ITYPE. EQ.l)  THEN 
DO  10  I  =  1,  IPARA 

AHI(I)  =  ANOM(I)*AKLO(I) 

ALO(I)  =  ANOM (I)  -  (AHI(I)  -  ANOM(I)) 

CONTINUE 

ELSE  IF  (ITYPE.EQ.2)  THEN 
DO  20  I  =  1,  IPARA 

ALO(I)  =  ANOM(I)/AKLO(I) 

AHI(I)  =  ANOM(I) *AKLO(I) 

CONTINUE 

ELSE  IF  (ITYPE. Eq. 3)  THEN 
DO  30  I  =  1,  IPARA 

ALO(I)  =  ANOM(I)/AKLC(I) 

AHI(I)  =  ANOM(I) *AKHI (I) 

CONTINUE 


ELSE  IF  (ITYPE.EQ .4)  THEN 
DO  40  I  *  1,  IPARA 

ALO(I)  =  ANOM(I)/AKLO(I) 

AHI(I)  =  ANOM(I) *AKHI (I) 

40  CONTINUE 
END  IF 
C 

C  ALGORITHM  EVALUATION  FOR  CASE  OF  NORMALLY  DISTRIBUTED  UNCERTAINTY. 

C 

IF  (ITYPE.EQ. 1)  THEN 
C 

C  EVALUATE  THE  ALGORITHM  THE  USER  REQUESTED  (ISIZE)  NUMBER  OF  TIMES. 

C 

DO  60  I  *  1,  ISIZE 
C 

C  RANDOMIZE  THE  PARAMETER  VALUES. 

C 

DO  50  J  =  1,  IPARA 
CALL  FNRN (NUMNR , RN ) 

A(J)  =  ANOM(J)  +  RN*(AHI( J)  -  AN0M(J))/1.96 
50  CONTINUE 

C 

C  EVALUATE  THE  ALGORITHM  USING  THE  RANDOMIZED  PARAMETER  VALUES. 

C 

CALL  TFCANOM.A.XVAL.YVAL) 

C 

C  APPLY  ALGORITHMIC  UNCERTAINTY  IF  REQUIRED  (I ALGO  =  1) 

C  AND  STORE  THE  RESULT  OF  THE  ALGORITHM  EVALUATION. 

C 

IF  (IALGO.NE.l)  THEN 
RESULT ( I )  =  YVAL 
ELSE 

CALL  FNRN (NUMNR ,RN) 

AHI(IPARAM)  =  YVAL*AKLO(IPARAM) 

RESULT (I)  =  YVAL  +  RN* (AHI (IPARAM)  -  YVAL)/1.96 
END  IF 
60  CONTINUE 
RETURN 
C 

C  ALGORITHM  EVALUATION  FOR  CASE  OF  LOGNORMALLY  DISTRIBUTED  UNCERTAINTY. 

C 

ELSE  IF  (ITYPE.EQ. 2)  THEN 
C 

C  EVALUATE  THE  ALGORITHM  THE  USER  REQUESTED  (ISIZE)  NUMBER  OF  TIMES. 

C 


DO  120  I  =  1,  ISIZE 
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C  RANDOMIZE  THE  PARAMETER  VALUES. 

C 

DO  110  J  -  1,  IPARA 
CALL  FNRN (NUMNR,RN) 

A(J)  =  ANOM(J) *AKLO ( J) ** (RN/ 1 . 96) 

110  CONTINUE 

C 

C  EVALUATE  THE  ALGORITHM  USING  THE  RANDOMIZED  PARAMETER  VALUES. 

C 

CALL  TF ( ANOM , A , XVAL , YVAL) 

C 

C  APPLY  ALGORITHMIC  UNCERTAINTY  IF  REQUIRED  (I ALGO  =  1) 

C  AND  STORE  THE  RESULT  OF  THE  ALGORITHM  EVALUATION. 

C 

IF  (IALGO.NE.l)  THEN 
RESULT (I)  *  YVAL 
ELSE 

CALL  FNRN (NUMNR , RN ) 

RESULT (I)  =  YVAL* AKLO(IPARAM)**(RN/ 1.96) 

END  IF 
120  CONTINUE 
RETURN 
C 

C  ALGORITHM  EVALUATION  FOR  CASE  OF  BETA  DISTRIBUTED  UNCERTAINTY. 

C 

ELSE  IF  (ITYPE. EQ .3)  THEN 
C 

C  EVALUATE  THE  ALGORITHM  THE  USER  REQUESTED  (ISIZE)  NUMBER  OF  TIMES. 

C 

DO  220  I  =  1,  ISIZE 
C 

C  RANDOMIZE  THE  PARAMETER  VALUES . 

C 

DO  210  J  *  1,  IPARA 
CALL  UPR1 (NUMUR.RN) 

LL  =  IBD(J) 

FL1  =  RN*BDMAX+1.0 
L  =  FL1 
LP1  =  L+l 
FL2  =  L 

FL1  =  FL1  -  FL2 
FL2  =  1.0  -  FL1 

BDV  =  FL1*BD(LP1 ,LL)  +  FL2*BD(L , LL) 

A(J)  =  ALO(J)  +  BDV* (AHI ( J) -ALO( J) ) 

210  CONTINUE 
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c 

C  EVALUATE  THE  ALGORITHM  USING  THE  RANDOMIZED  PARAMETER  VALUES . 

C 

CALL  TF(ANOM,A,XVAL, YVAL) 

C 

C  APPLY  ALGORITHMIC  UNCERTAINTY  IF  REQUIRED  (I ALGO  =  1) 

C  AND  STORE  THE  RESULT  OF  THE  ALGORITHM  EVALUATION. 

C 

IF  (IALGO.NE.l)  THEN 
RESULT(I)  =  YVAL 
ELSE 

CALL  UPR1 (NUMUR , RN ) 

ALO(IPARAM)  =  YVAL/AKLO( IPARAM) 

AHI (IPARAM)  =  YVAL+AKHI (IPARAM) 

LL  =  IBD (IPARAM) 

FL1  =  RN*BDMAX+1 . 0 
L  =  FL1 
LP1  =  L+l 
FL2  =  L 

FL1  =  FL1  -  FL2 
FL2  *  1.0  -  FL1 

BDV  *  FL1*BD(LP1,LL)  +  FL2*BD(L,LL) 

RESULT (I)  =  ALO (IPARAM)  +  BDV* (AHI (IPARAM)  -  ALO (IPARAM)) 
END  IF 
220  CONTINUE 
RETURN 
C 

C  ALGORITHM  EVALUATION  FOR  CASE  OF  UNIFORMLY  DISTRIBUTED  UNCERTAINTY. 

C 

ELSE  IF  (ITYPE. EQ .4)  THEN 
C 

C  EVALUATE  THE  ALGORITHM  THE  USER  REQUESTED  (ISIZE)  NUMBER  OF  TIMES. 

C 

DO  320  I  =  1,  ISIZE 
C 

C  RANDOMIZE  THE  PARAMETER  VALUES. 

C 

DO  310  J  =  1,  IPARA 
CALL  UPR1 (NUMUR, RN) 

IF  (RN.LE.0.5)  THEN 

A(J)  =  ALO(J)  +  2 . 0*RN* (ANOM( J)  -  ALO(J)) 

ELSE 

A( J)  =  ANOM(J)  +  2 . 0* (RN  -  0.5)*(AHI(J)  -  ANOM(J)) 

END  IF 

310  CONTINUE 


C 


G  EVALUATE  THE  ALGORITHM  USING  THE  RANDOMIZED  PARAMETER  VALUES. 

C 

CALL  TF (ANOM , A ,XVAL , YVAL) 

C 

C  APPLY  ALGORITHMIC  UNCERTAINTY  IF  REQUIRED  (I ALGO  =  1) 

C  AND  STORE  THE  RESULT  OF  THE  ALGORITHM  EVALUATION. 

C 

IF  (IALGO.NE.l)  THEN 
RESULT (I)  =  YVAL 
ELSE 

CALL  UPR1 (NUMUR.RN) 

ALO(IPARAM)  *  YVAL/AKLO(IPARAM) 

AHI (IPARAM)  =  YVAL*AKHI( IPARAM) 

IF  (RN.LE.O .5)  THEN 

RESULT (I)  =  ALO (IPARAM)  +  2.0*RN*(YVAL  -  ALO (IPARAM)) 
ELSE 

RESULT (I)  =  YVAL  +  2.0*(RN  -  0 . 5)* (AHI (IPARAM)  -  YVAL) 
ENDIF 
END  IF 
320  CONTINUE 
RETURN 
ENDIF 
RETURN 
END 
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SUBROUTINE  sort(x,inum) 


C  SUBROUTINE  SORT  SORTS  THE  INUM  VALUES  IN  THE  ARRAY  X  FROM  LOWEST 

C  TO  HIGHEST  VALUE.  IN  THE  PRESENT  CONTEXT,  THE  SUBROUTINE  IS  USED 

C  TO  SORT  THE  ARRAY  OF  COMPUTED  ALGORITHM  RESULTS  FROM  LOWEST  TO 

C  HIGHEST  VALUE. 

C 

INTEGER+4  INUM 
REAL*4  X(*) 

NEND  =  INUM  -  1 
10  IDONE  =  0 

DO  20  I  *  1,  NEND 

IF  (X(I) .GT.XCl+l))  THEN 
XTEM  =  X(I) 

X(I)  =  X(I+1) 

X(I+1)  =  XTEM 
IDONE  =  1 
END  IF 
20  CONTINUE 

IF  (IDONE. Eq.O)  GO  TO  30 
NEND  =  NEND  -  1 
GO  TO  10 
30  CONTINUE 
RETURN 
END 
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SUBROUTINE  STATISTIC(X,  XINT,  XOUT,  XXLN,  XXII,  I FREQ) 


G 

C  SUBROUTINE  STATISTIC  CALCULATES  VARIOUS  STATISTICS  FOR  THE  ORDERED 

C  (FROM  LOWEST  TO  HIGHEST  VALUE)  ALGORITHM  RESULTS  DISTRIBUTION  WHICH 

C  IS  PASSED  TO  THE  SUBROUTINE  IN  ARRAY  X.  FOUR  BASIC  MANIPULATIONS  OF 

C  THE  ALGORITHM  RESULTS  DISTRIBUTION  ARE  PERFORMED  BY  THE  SUBROUTINE: 

C 

C  1.  GENERATE  HISTOGRAM  DATA. 

C 

C  THE  FULL  RANGE  (FROM  X(l)  TO  X(ISIZE)  WHERE  ISIZE  IS  THE  NUMBER 

C  OF  COMPUTED  ALGORITHM  RESULTS)  OF  THE  RESULTS  DISTRIBUTION  IS 

C  DIVIDED  INTO  NINTV  EQUAL  LENGTH  INTERVALS,  THE  BOUNDARY  VALUES 

C  OF  WHICH  ARE  CONTAINED  IN  THE  ARRAY  XINT  (NINTV+1  ENTRIES) .  THE 

C  NUMBER  OF  RESULTS  DISTRIBUTION  VALUES  OCCURRING  IN  EACH  RANGE 

C  INTERVAL  IS  THEN  DETERMINED  AND  STORED  IN  THE  INTEGER  ARRAY  IFREQ . 

C  THE  RESULT  IS  A  HISTOGRAM  FOR  THE  RESULTS  DISTRIBUTION. 

C 

C  2 .  DETERMINE  DATA  MEAN  (OR  AVERAGE) . 

C 

C  FOR  NORMALLY,  LOGNORMALLY,  AND  BETA  DISTRIBUTED  ALGORITHM  RESULTS, 

C  THE  MIDPOINT  OF  THE  I-TH  RANGE  SUBINTERVAL  IS  MULTIPLIED  BY  THE 

C  NUMBER  OF  OCCURRENCES  OF  THE  ALGORITHM  VALUE  IN  THE  SUBINTERVAL; 

C 

0 . 5*(XINT(I)  -  XINT(I+1))*IFREQ(I) 

C  AND  THE  RESULT  IS  SUMMED  OVER  ALL  SUBINTERVALS .  THE  EVALUATED 

C  SUM  IS  THEN  DIVIDED  BY  THE  TOTAL  NUMBER  OF  ALGORITHM  RESULTS  (ISIZE) 

C  TO  GENERATE  THE  APPROXIMATE  MEAN  OF  THE  RESULTS  DISTRIBUTION. 

C 

C  FOR  UNIFORMLY  DISTRIBUTED  ALGORITHM  RESULTS,  THE  FRACTIONAL  POSITION 

C  IN  THE  INTERVAL  CORRESPONDING  TO  THE  ISIZE/2  RESULT  IN  THE  ORDERED 

C  RESULTS  DISTRIBUTION  IS  DETERMINED  AND  REPORTED  AS  THE  MEAN  (AMEAN) 

C  OF  THE  DISTRIBUTION. 

C 

C  3.  ESTIMATE  VARIANCE  AND  STANDARD  DEVIATION: 

C 

C  THE  COMPUTED  MEAN  OF  THE  DISTRIBUTION  IS  SUBTRACTED  FROM  THE  MIDPOIN 

C  VALUE  (XAVE)  OF  EACH  SUBINTERVAL  TO  DEFINE  VARIABLE  XDIF: 

C 

C  XDIF  »  0 . 5*(XINT(I+1)  +  XINT(I))  -  AMEAN  =  XAVE  -  AMEAN 

C 

C  THE  VALUE  OF  XDIF  FOR  EACH  SUBINTERVAL  IS  NEXT  SQUARED,  MULTIPLIED 

C  BY  THE  NUMBER  OF  ALGORITHM  RESULTS  OCCURRING  IN  THE  SUBINTERVAL,  AND 

C  THE  RESULT  FOR  EACH  SUBINTERVAL  IS  SUMMED  OVER  ALL  THE  INTERVALS. 

C  THE  COMPUTED  SUM  IS  THEN  DIVIDED  BY  THE  TOTAL  NUMBER  OF  VALUES  (ISIZE) 
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C  IN  THE  RESULTS  DISTRIBUTION,  AND  THE  RESULT  (VAR)  IS  INTERPRETED  AS 

C  AN  APPROXIMATE  VALUE  FOR  THE  VARIANCE  OF  THE  ALGORITHM  RESULTS 

C  DISTRIBUTION.  HIE  SQUARE  ROOT  OF  VAR  THEN  PROVIDES  THE  ESTIMATE  FOR 

C  THE  STANDARD  DEVIATION  (SIGMA)  OF  THE  RESULTS  DISTRIBUTION. 

C 

C  4.  ESTIMATE  THE  DISTRIBUTION  K-FACTORS: 

C 

C  NORMALLY  DISTRIBUTED  RESULTS: 

C 

C  K  =  1.0  +  1 . 96*SIGMA/AMEAN 

LOGNORMALLY  DISTRIBUTED  RESULTS: 

C 

C  K  =  EXP(1 . 96*SIGMA/AMEAN) 

C 

C  BETA  AND  UNIFORMLY  DISTRIBUTED  RESULTS: 

C 

C  KLO  =  AMEAN/X(1) 

C  KHI  =  X(ISIZE)/AMEAN 

C  KAVE  =  0.5* (KLO  +  KHI) 

C 

C  THE  SUBROUTINE  THEN  INQUIRES  WHETHER  OR  NOT  PLOT  FILES  FOR  THE 

C  ALGORITHM  RESULTS  DISTRIBUTION  AND  FITS  TO  THE  RESULTS  ARE  TO  BE 

C  GENERATED.  IF  GENERATION  OF  PLOT  FILES  IS  REQUESTED  BY  THE  USER, 

C  THE  SUBROUTINE  CALLS  SUBROUTINE  HISTO. 

C 

COMMON  /BLKOl/  IPARAM,  I TYPE,  ISIZE,  I SEED,  IALGO ,  IXVAL,  ITABL , 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

REAL*4  X(*),  XINT(*) ,  XOUT(*) ,  XXLN(*) ,  XXII(*) 

INTEGER*4  IFREQ(*) 

CHARACTER* 1  ANS2 
C 

C  COMPUTE  INCREMENT  FOR  RESULTS  SUBINTERVALS. 

C 

DEL  *  (X(ISIZE)  -  X(l) ) /NINTV 
C 

C  GENERATE  RESULTS  SUBINTERVALS. 

C 

XINT(l)  *  X(l) 

XINT(NINTV+1)  =  X(ISIZE) 

DO  10  I  -  2,  NINTV 

XINT(I)  =  XINT(I-l)  +  DEL 
10  CONTINUE 
C 

C  DETERMINE  THE  NUMBER  OF  ALGORITHM  RESULTS  OCCURRING  IN  EACH 

C  RANGE  SUBINTERVAL  AND  STORE  IN  ARRAY  I FREQ . 
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DO  20  I  *  1,  NINTV 
IFREQ(I)  *  0 
20  CONTINUE 

DO  40  I  *  1,  ISIZE 
RES  =  X(I) 

DO  30  J  *  1,  NINTV 

IF  (RES.GT.XINT(J+1))  THEN 
GO  TO  30 

ELSE  IF  ((RES.GE.XINT(J)) .AND. (RES.LT.XINT(J+1)))  THEN 
IFREq(J)  =  IFREQ(J)  +  1 
GO  TO  40 
ENDIF 
30  CONTINUE 
40  CONTINUE 

IFREq (NINTV)  =  IFREq (NINTV)  +  1 
C 

C  COMPUTE  THE  DISTRIBUTION  MEAN  FOR  NORMALLY,  LQGNORMALLY,  AND 
C  BETA  DISTRIBUTED  ALGORITHM  RESULTS  DISTRIBUTIONS. 

C 

IF  ( (ITYPE.Eq .1) .OR. (ITYPE.Eq.2) . OR. (ITYPE. Eq.3) )  THEN 
AMEAN  =0.0 
DO  50  I  -  1,  NINTV 

AMEAN  =  AMEAN  +  0 . 5* (XINT(I+l)  +  XINT(I))*IFREq(I) 

50  CONTINUE 

AMEAN  =  AMEAN/ ISIZE 
GO  TO  70 
C 

C  COMPUTE  THE  DISTRIBUTION  MEAN  FOR  UNIFORMLY  DISTRIBUTED 
C  ALGORITHM  RESULTS  DISTRIBUTIONS. 

C 

ELSE  IF  (ITYPE.Eq. 4)  THEN 
IRESULT  =  0 
DO  60  I  =  1,  NINTV 

IDIF  *  ISIZE/2  -  IRESULT  -  TFREq(I) 

IF  (IDIF.GT.O)  THEN 

IRESULT  =  IRESULT  +  IFREq (I) 

ELSE  IF  (IDIF.LE.O)  THEN 

FRAC  =  (ISIZE/2  -  IRESULT) /IFREq (I) 

AMEAN  =  XINT(I)  +  FRAC*(XINT(I+1)  -  XINT(I)) 

GO  TO  70 
ENDIF 
60  CONTINUE 
ENDIF 

70  VAR  =  0.0 
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c 

c 
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c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


COMPUTE  THE  DISTRIBUTION  VARIANCE  FOR  ALL  ALGORITHM  RESULTS 
DISTRIBUTION  TYPES. 

DO  80  I  =  1,  NINTV 

XAVE  =  0.5*(XINT(I+1)  +  XINT(I)) 

XDIF  -  XAVE  -  AMEAN 

VAR  =  VAR  +  XDIF  *XDIF  * IFREQ ( I ) 

80  CONTINUE 

VAR  *  VAR/ISIZE 

COMPUTE  THE  STANDARD  DEVIATION  (SQRT(VAR))  OF  THE  DISTRIBUTION. 

SIGMA  =  SQRT(VAR) 

OUTPUT  THE  COMPUTED  MEAN,  VARIANCE,  AND  STANDARD  DEVIATION  FOR 
THE  ALGORITHM  RESULTS  DISTRIBUTION. 

WRITE (NPRT, 1000)  AMEAN,  VAR,  SIGMA,  ISIZE 

WRITE (NOUT, 1000)  AMEAN,  VAR,  SIGMA,  ISIZE 

COMPUTE  AND  OUTPUT  THE  K-FACTOR  FOR  NORMALLY  DISTRIBUTED 
ALGORITHM  RESULTS. 

IF  (ITYPE.EQ. 1)  THEN 

AKEST  -  1.0  +  1 . 96*SIGMA/ AMEAN 
WRITE(NPRT, 1100)  AKEST 
WRITE(NOUT, 1100)  AKEST 

COMPUTE  AND  OUTPUT  THE  K-FACTOR  FOR  LOGNORMALLY  DISTRIBUTED 
ALGORITHM  RESULTS. 

ELSE  IF  (ITYPE.EQ. 2)  THEN 

AKEST  =  EXP ( 1 . 96*SIGMA/AMEAN ) 

WRITE(NPRT, 1100)  AKEST 
WRITE(NQUT, 1100)  AKEST 

COMPUTE  AND  OUTPUT  THE  K-FACTORS  (KLO,  KHI ,  KAVE)  FOR  BETA  AND 
UNIFORMLY  DISTRIBUTED  ALGORITHM  RESULTS. 

ELSE  IF  ((ITYPE.EQ. 3) .OR. (ITYPE.EQ. 4))  THEN 
AKEST1  =  AMEAN/X(1) 

AKEST2  =  X(ISIZE) /AMEAN 
AKAVE  =  0 . 5*(AKEST1  +  AKEST2) 

WRITE(NPRT, 1200)  AKEST1 ,  AKEST2 ,  AKAVE 
WRITE(NOUT, 1200)  AKEST 1 ,  AKEST2 ,  AKAVE 

END  IF 
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C  DETERMINE  IF  A  RESULTS  DISTRIBUTION  PLOTFILE  AND  A  FIT  TO  THE 
C  ALGORITHM  RESULTS  DISTRIBUTION  ARE  TO  BE  GENERATED.  IF  SO, 

C  CALL  SUBROUTINE  HISTO,  ELSE  EXIT  SUBROUTINE. 

C 

WRITE (NOUT, 1300) 

READ(NINP, 1031)  ANS2 
1031  FORMAT (A 1) 

IF  ((ANS2.EQ.  ’Y’).0R.(ANS2.EQ.  'y'))  THEN 

CALL  HISTO (XINT,  I FREQ ,  AMEAN,  SIGMA,  XOUT,  XXLN,  XXII) 

END  IF 
RETURN 

1000  FORMAT(/,’  MEAN  OF  DISTRIBUTION  =  \E12.5,/, 

1  '  VARIANCE  OF  DISTRIBUTION  =  \E12.5,/, 

2  '  STANDARD  DEVIATION  [SQRT (VARIANCE)]  =  ’,E12.5,/, 

3  '  TOTAL  NUMBER  OF  SAMPLES  =  '.6X.I6) 

1100  FORMAT (/ , '  K-FACTOR  ESTIMATE:  K  =  \F7.4,/) 

1200  FORMAT(/,’  K-FACTOR  ESTIMATES:’,/,/, 

1  ’  K-LOW  =  ’ ,F7.4,/, 

2  ’  K-HIGH  =  ’ ,F7.4,/, 

3  ’  K-AVERAGE  =  ’,F7.4,/) 

1300  FORMAT (/,'  DO  YOU  WANT  PLOT  FILES  OF  THE  RESULTS’,/, 

1  ’  DISTRIBUTION  AND  FIT  TO  BE  GENERATED?  [Y/N]',/) 

END 


SUBROUTINE  HISTO(XX,  II,  AMEAN,  SIGMA,  XOUT,  XXLN,  XXII) 

C 

C  SUBROUTINE  HISTO  GENERATES  (X,Y)  DATA  FILES  WHICH  CAN  BE  USED 

C  TO  GENERATE  PLOTS  OF  THE  ALGORITHM  RESULTS  DISTRIBUTION, 

C  COMPUTED  FITS  TO  THE  DISTRIBUTION  RESULTS,  AND  THE  CUMULATIVE 

C  PROBABILITY  DISTRIBUTION  FOR  THE  RESULTS.  THE  (X,Y)  PLOT  DATA 
C  ARE  WRITTEN  TO  DEFAULT  UNIT  3. 

C 

C  THE  BASIC  DATA  DELIVERED  TO  THE  SUBROUTINE  ARE  AS  FOLLOWS: 

C 

C  XX  ARRAY  OF  DISTRIBUTION  SUBINTERVAL  BOUNDARIES. 

C  II  INTEGER  ARRAY  OF  NUMBER  OF  ALGORITHM  RESULTS 

C  OCCURRING  IN  EACH  SUBINTERVAL. 

C  AMEAN  -  COMPUTED  MEAN  (OR  AVERAGE)  OF  THE  DISTRIBUTION. 

C  SIGMA  -  COMPUTED  STANDARD  DEVIATION  OF  THE  DISTRIBUTION. 

C 

C  THE  SUBROUTINE  PERFORMS  THE  FOLLOWING  OPERATIONS: 

C 

C  1.  OUTPUT  HISTOGRAM  DATA. 

C 

C  THE  MIDPOINT  VALUE  XOUT  OF  THE  I-TH  SUBINTERVAL  IS  COMPUTED 

C 

C  XOUT(I)  *  0.5*(XX(I)  +  XXCI+D) 

C 

C  AND  THE  DATA  PAIRS  (XOUT(I) ,  II(I))  FOR  EACH  SUBINTERVAL  ARE 

C  OUTPUT  TO  THE  PLOT  FILE. 

C 

C  2.  CURVE  FIT  TO  HISTOGRAM  DATA. 

C 

C  A.  NORMALLY  DISTRIBUTED  ALGORITHM  RESULTS: 

C 

C  THE  SUBINTERVAL  CONTAINING  THE  LARGEST  NUMBER  OF  COMPUTED 

C  ALGORITHM  RESULTS  IS  DETERMINED  AND  INTEGER  VARIABLE  UMAX 

C  IS  SET  EqUAL  TO  THE  NUMBER  OF  ALGORITHM  RESULTS  WHICH  OCCURRED 

C  IN  THAT  SUBINTERVAL.  NEXT,  FOR  EACH  SUBINTERVAL  THE  DIFFERENCE 

C  BETWEEN  THE  MIDPOINT  OF  THE  SUBINTERVAL  AND  THE  DISTRIBUTION 

C  MEAN  IS  EVALUATED  AND  DIVIDED  BY  THE  STANDARD  DEVIATION  TO 

C  PRODUCE  VARIABLE  TEM: 

C 

C  TEM  =  (XOUT(I)  -  AMEAN) /SIGMA 

C 

C  THEN  THE  NORMAL  DISTRIBUTION  FIT  FOR  THE  I-TH  INTERVAL  IS 

C  COMPUTED; 

C 

C  TEM  =  IIMAX*EXP(-0.5*TEM*TEM) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


AND  THE  DATA  PAIR  (XOUT(I),  TEM)  IS  OUTPUT  TO  THE  PLOT  FILE. 

B.  LOGNORMALLY  DISTRIBUTED  ALGORITHM  RESULTS: 

THE  RATIO  SSIGMA  =  SIGMA/AMEAN  IS  COMPUTED  AND  THE  NATURAL 
LOGARITHM  OF  AMEAN  IS  SET  EQUAL  TO  VARIABLE  AMLN.  NEXT  THE 
NATURAL  LOGARITHM  OF  EACH  SUBINTERVAL  BOUNDARY  IS  COMPUTED 
AND  STORED  IN  ARRAY  XXLN.  THEN  THE  SUBINTERVAL  CONTAINING 
THE  GREATEST  NUMBER  OF  ALGORITHM  RESULTS  IS  DETERMINED  AND 
INTEGER  VARIABLE  UMAX  IS  SET  EQUAL  TO  THE  NUMBER  OF  ALGORITHM 
RESULTS  WHICH  OCCURRED  IN  THAT  SUBINTERVAL.  NEXT  A  NORMAL 
CURVE  VALUE  (TEM1  AND  TEM2)  IS  COMPUTED  AT  THE  BOUNDARY  OF 
EACH  SUB INTERVAL; 

TEM1  =  (XXLN(I)  -  AMLN) /SSIGMA 
TEM1  =  EXP(-0.5*TEM1*TEM1) 

TEM2  =  (XXLN(I+1)  -  AMLN) /SSIGMA 
TEM2  =  EXP(-0.5*TEM2*TEM2) 

AND  THE  AVERAGE  OF  THE  TWO  RESULTS  (TEM)  IS  DETERMINED. 

MULTIPLIED  BY  THE  WIDTH  OF  THE  SUBINTERVAL  (XXLNQ+i)  -  XXLN (I)), 
AND  STORED  IN  ARRAY  XXII.  THE  RESULTING  LOGNORMAL  DISTRIBUTION 
IS  THEN  SCALED  BY  UMAX  AND  OUTPUT  TO  THE  PLOT  FILE. 

C.  BETA  DISTRIBUTED  ALGORITHM  RESULTS: 

THE  INTEGER  VARIABLE  UMAX  IS  OBTAINED  AS  ABOVE  FOR  THE  NORMALLY 
DISTRIBUTED  CASE.  NEXT  THE  COMPUTED  MEAN  AND  STANDARD  DEVIATION 
OF  THE  RESULTS  DISTRIBUTION  ARE  SUPPLIED  TO  SUBROUTINE  BETAFIT 
IN  WHICH  AN  APPROPRIATE  TABULATED  OR  CALCULATED  BETA  DISTRIBUTION. 
CORRESPONDING  TO  THE  COMPUTED  MEAN  AND  STANDARD  DEVIATION,  IS 
GENERATED.  THIS  GENERATED  CURVE  IS  THEN  SCALED  USING  THE  VALUE 
OF  UMAX  AND  THE  RESULTING  DISTRIBUTION  IS  OUTPUT  TO  THE  PLOT  FILE 

D.  UNIFORMLY  DISTRIBUTED  ALGORITHM  RESULTS: 

THE  COMPUTED  MEAN  OF  THE  ALGORITHM  RESULTS  DISTRIBUTION  IS 
COMPARED  TO  THE  FULL  RANGE  (XX(1)  TO  XX(NINTV+l))  OF  THE 
DISTRIBUTION  TO  DETERMINE  ITS  RELATIVE  LOCATION  IN  THE  FULL 
RANGE.  A  UNIFORM  DISTRIBUTION  WITH  507.  OF  ITS  VALUES  BETWEEN 
XX (1)  AND  AMEAN  AND  507.  BETWEEN  AMEAN  AND  XX(NINTV+l)  IS  THEN 
COMPUTED  AND  OUTPUT  TO  THE  PLOT  FILE. 

COMMON  /BLK01/  IPARAM,  I TYPE,  ISIZE,  I SEED,  I ALGO ,  IXVAL,  ITABL , 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 
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REAL*4  XX(*),  XOUT (*)  ,  XXLN(*) ,  XXIIO) 

INTEGER*4  II (*) 

C 

C  EVALUATE  THE  SUBINTERVAL  MIDPOINTS  AND  OUTPUT  THE  HISTOGRAM 

C  DATA  FOR  THE  COMPUTED  RESULTS  DISTRIBUTION  TO  THE  PLOT  FILE. 

r\ 

DO  10  I  3  1,  NINTV 

XOUT(I)  *  0 . 5*(XX(I)  +  XX(I+D) 

WRITE (NPLT, 1000)  XOUT (I),  II (I) 

10  CONTINUE 
C 

C  COMPUTE  A  FIT  TO  NORMALLY  DISTRIBUTED  ALGORITHM  RESULTS. 

C 

IF  (ITYPE.EQ. 1)  THEN 
UMAX  =  0 

DO  100  I  =  1,  NINTV 

IF  (II (I)  .GT. UMAX)  UMAX  =  II (I) 

100  CONTINUE 

DO  110  I  =  1,  NINTV 

TEM  =  (XQUT(I)  -  AMEAN) /SIGMA 
TEM  «  IIMAX*EXP(-0 . 5*TEM*TEM) 

WRITE(NPLT.llOO)  XOUT(I) ,  TEM 
110  CONTINUE 
C 

C  COMPUTE  A  FIT  TO  LOGNORMALLY  DISTRIBUTED  ALGORITHM  RESULTS. 

C 

ELSE  IF  (ITYPE.EQ. 2)  THEN 
SSIGMA  =  SIGMA/AMEAN 
AMLN  =  ALOG (AMEAN) 

DO  200  I  =  1,  NINTV+1 
XXLN(I)  *  ALOG (XX (I) ) 

200  CONTINUE 
UMAX  -  0 

DO  210  I  =  1,  NINTV 

IF  ((XX(I)  .LT. AMEAN)  .AND.  (XX(I+1)  .GE. AMEAN))  UMAX  =  II(I) 
210  CONTINUE 

DO  220  I  =  1,  NINTV 

TEM1  =  (XXLN(I)  -  AMLN) /SSIGMA 
TEM1  =  EXP(-0.5*TEM1*TEM1) 

TEM2  =  (XXLN(I+1)  -  AMLN) /SSIGMA 
TEM2  =  EXP(-0.5*TEM2*TEM2) 

TEM  *  0 . 5*(TEM1  +  TEM2) 

XXII(I)  =  TEM*(XXLN(I+1)  -  XXLN(I)) 

220  CONTINUE 

DO  230  I  -  1,  NINTV 

IF  ( (XXLN(I) . LT. AMLN) .AND. (XXLN(I+1) . GE.AMLN)) 
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1  SCALE  =  IIMAX/XXlKD 

230  CONTINUE 

DO  240  I  ■  1,  NINTV 
TEM  =  SCALE*XXII (I) 

WRITE(NPLT.llOO)  XOUT(I) ,  TEM 
240  CONTINUE 
C 

C  COMPUTE  A  FIT  TO  BETA  DISTRIBUTED  ALGORITHM  RESULTS. 

C 

ELSE  IF  (ITYPE.EQ .3)  THEN 
UMAX  =  0 

DO  300  I  *  1,  NINTV 

IF  ((XX(I) .LT.AMEAN) .AND. (XX(I+1) .GE.AMEAN))  UMAX  =  I 
300  CONTINUE 
C 

C  CALL  BETAFIT  TO  EVALUATE  THE  ACTUAL  BETA  DISTRIBUTION  FIT. 

C 

CALL  BETAFIT (AMEAN,  SIGMA,  XX,  XXLN) 

SCALE  =0.0 
NSCALE  =  5 

DO  310  I  =  UMAX  -  NSCALE,  UMAX  +  NSCALE 

SCALE  =  SCALE  +  2 . 0*11 (I)/(XXLN(I)  +  XXLN(I+1)) 

310  CONTINUE 

SCALE  =  SCALE/ (2*NSCALE  +  1) 

DO  320  I  =  1,  NINTV 

TEM  =  0.5*(XXLN(I)  +  XXLN(I+1) ) *SCALE 
WRITE (NPLT, 1100)  XOUT(I),  TEM 
320  CONTINUE 
C 

C  COMPUTE  A  FIT  TO  UNIFORMLY  DISTRIBUTED  ALGORITHM  RESULTS. 

C 

ELSE  IF  (ITYPE.EQ .4)  THEN 

SCALE  =  NINTV* (AMEAN  -  XX(1))/(XX(NINTV+1)  -  XX(l)) 

DO  400  I  =  1,  NINTV 

IF  (XOUT(I) .LE.AMEAN)  THEN 
TEM  =  0 . 5*ISIZE/SCALE 
ELSE 

TEM  =  0 . 5*ISIZE/ (NINTV  -  SCALE) 

END  IF 

WRITE (NPLT, 1100)  XOUT(I) ,  TEM 
400  CONTINUE 
END  IF 
C 

C  COMPUTE  THE  CUMULATIVE  PROBABILITY  DISTRIBUTION  ASSOCIATED 

C  WITH  THE  COMPUTED  ALGORITHM  RESULTS  DISTRIBUTION. 

C 
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WRITE (NPLT, 1100)  XX(1),  0.0 

AREA  =0.0 

DO  410  I  =  1,  NINTV 

AREA  =  AREA  +  (XX(I+1)  -  XX(I))*II(I) 
410  CONTINUE 

TOTAREA  =  AREA 

AREA  =0.0 

DO  420  I  =  1,  NINTV 

AREA  =  AREA  +  (XX(I+1)  -  XX(I))*II(I) 
WRITE (NPLT, 1100)  XX(I+1),  AREA/TOTAREA 
420  CONTINUE 
RETURN 

1000  FORMAT (E12. 5, ' ,  ',15) 

1100  F0RMAT(E12.5, ’ ,  \E12.5) 

END 
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C  SIGMA  =  SIGMA/ (XX (NINTV+1)  -  X(I)) 

C  TEM  (AMEDIAN  -  AMEDIAN*AMEDIAN)/ (SIGMA*SIGMA)  -  1. 

C  ALPHA  =  AMEDIAN+TEM 

C  BETA  =  (1.0  -  AMEDIAN) *TEM 

C 

C  ONCE  THE  VALUES  FOR  ALPHA  AND  BETA  ARE  AVAILABLE,  FUNCTION  ROUTINE  FX 

C  IS  CALLED  TO  COMPUTE  THE  (UNNORMALIZED)  BETA  DISTRIBUTION,  THE  AREA 

C  UNDER  THE  DISTRIBUTION  IS  COMPUTED,  AND  THEN  THE  DISTRIBUTION  IS 

C  NORMALIZED.  THE  RESULTING  DISTRIBUTION  IS  RETURNED  TO  THE  CALLING 

C  SUBROUTINE  IN  ARRAY  Y. 

C 

COMMON  /BLK01/  IPARAM,  I TYPE,  ISIZE,  I SEED,  I ALGO ,  IXVAL,  ITABL, 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

REAL+4  AMED(81) ,  ALPH(81) ,  BET(8l) 

REAL*4  XX(*),  Y(*) 

DATA  AMED/O. 167, 0.172, 0.179, 0.185, 0.192, 0.200, 0.208, 0.217, 0.227, 

1  0.238,0.250,0.256,0.263,0.270,0.278,0.286,0.294,0.303, 

2  0.313,0.323,0.333,0.339,0.345,0.351,0.357,0.364,0.370, 

3  0.377,0.385,0.392,0.400,0.408,0.417,0.426,0.435,0.444, 

4  0.455,0.465,0.476,0.488,0.500,0.512,0.524,0.535,0.545, 

5  0.556,0.565,0.574,0.583,0.592,0.600,0.608,0.615,0.623, 

6  0.630,0.636,0.643,0.649,0.655,0.661,0.667,0.677,0.687, 

7  0.697,0.706,0.714,0.722,0.730,0.737,0.744,0.750,0.762, 

8  0.773,0.783,0.792,0.800,0.808,0.815,0.821,0.828,0.833/ 


DATA 

ALPH/ 1.04, 

1.07, 

1.11, 

1.14, 

1.18, 

1.23, 

1.28, 

1.33, 

1.39, 

1 

1.46, 

1.54, 

1.58, 

1.63, 

1.68, 

1.73, 

1.79, 

1.85, 

1.92, 

2 

1.99, 

2.07, 

2.16, 

2.20, 

2.26, 

2.31, 

2.36, 

2.42, 

2.48, 

3 

2.55, 

2.62, 

2.69, 

2.77, 

2.86, 

2.95, 

3.05, 

3.15, 

3.26, 

4 

3.39, 

3.52, 

3.66, 

3.83, 

41*4.0/ 

DATA 

BET/41+4.0 

,  3.83, 

3.66 

,3.52 

,3.39 

> 

1 

3.26, 

3.15, 

3.05, 

2.95, 

2.86, 

2.77, 

2.69, 

2.62, 

2.55, 

2 

2.48, 

2.42, 

2.36, 

2.31, 

2.26, 

2.20, 

2.16, 

2.07, 

1.99, 

3 

1.92, 

1.85, 

1.79, 

1.73, 

1.68, 

1.63, 

1.58, 

1.54, 

1.46, 

4 

1.39, 

1.33, 

1.28, 

1.23, 

1.18, 

1.14, 

1.11, 

1.07, 

1.04/ 

C 

C  COMPUTE  THE  MEDIAN  OF  THE  ALGORITHM  RESULTS  DISTRIBUTION. 

C 

AMEDIAN  =  (AMEAN  -  XX(1))/(XX(NINTV  +  1)  -  XX(1)) 

C 

C  DETERMINE  ALPHA  AND  BETA  FOR  THE  TABULATED  DISTRIBUTION  CASE. 

C 

IF  (ITABL. EQ.l)  THEN 
TEMPI  =  0 . 5*(AMED(1)  +  AMED(2) ) 

TEMP 2  =  0 . 5*(AMED(80)  +  AMED(81) ) 

IF  (AMEDIAN. LE. TEMPI)  THEN 
IBET  =  1 
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ELSE  IF  ( AMEDI AN . GT . TEMP2 )  THEN 
IBET  =81 
ELSE 

DO  10  I  =  2,  80 

TEMPI  =  0 . 5*(AMED(I)  +  AMED(I-l)) 

TEMP 2  =  0 . 5*(AMED(I)  +  AMED(I+l)) 

IF  ((AMEDIAN. GT. TEMPI) .AND. (AMEDIAN.LE.TEMP2) )  IBET  =  I 
10  CONTINUE 
END  IF 

ALPHA  =  ALPH(IBET) 

BETA  =  BET (IBET) 

C 

C  DETERMINE  ALPHA  AND  BETA  FOR  THE  CALCULATED  DISTRIBUTION  CASE. 

C 

ELSE  IF  (ITABL.EQ.2)  THEN 

SIGMA  =  SIGMA/ (XX (NINTV  +  1)  -  XX(1)) 

TEM  =  (AMEDIAN  -  AMEDI AN* AMEDIAN) /(SIGMA*SIGMA)  -  1.0 
ALPHA  =  AMEDIAN*TEM 
BETA  =  (1.0  -  AMEDIAN) *TEM 
END  IF 
C 

C  INITIALIZE  THE  POINTS  IN  THE  UNIT  INTERVAL  AT  WHICH  THE  BETA 

C  DISTRIBUTION  WILL  BE  COMPUTED. 

C 

Y(l)  =  0.0 
Y (NINTV+1)  =1.0 
DO  20  I  =  2,  NINTV 

Y(I)  =  (XX(I)  -  XX(1))/(XX(NINTV  +  1)  -  XX(1)) 

20  CONTINUE 
C 

C  COMPUTE  THE  UNNORMALIZED  BETA  DISTRIBUTION  AND  ITS  AREA. 

C 

SUMNORM  =0.0 
DO  30  I  =  1,  NINTV 

FX1  =  FX (ALPHA,  BETA,  Y(I)) 

FX2  =  FX (ALPHA,  BETA,  Y(I+1)) 

SUMNORM  =  SUMNORM  +  0.5*(FX1  +  FX2)*(Y(I+1)  -  Y(I» 

30  CONTINUE 
C 

C  COMPUTE  THE  NORMALIZED  BETA  DISTRIBUTION. 

C 

DO  40  I  =  1,  NINTV+1 

FX1  =  FX( ALPHA,  BETA,  Y(I)) 

Y (I)  =  FX1/SUMN0RM 
40  CONTINUE 
RETURN 
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FUNCTION  FX( ALPHA,  BETA,  YY) 

C 

C  FUNCTION  FX  EVALUATES  THE  FUNCTION: 

C 

C  (ALPHA-1.)  (BETA-1.) 

C  YY  (l.-YY) 

C 

C  FOR  DELIVERED  VALUES  OF  ALPHA,  BETA,  AND  YY.  THE  ABOVE  FUNCTION 

C  IS  THE  INTEGRAND  IN  THE  INTEGRAL  FORMULA  FOR  THE  BETA 

C  PROBABILITY  DISTRIBUTION  WITH  PARAMETERS  ALPHA  AND  BETA. 

C 

IF  ((YY.GE.l.O) .OR. (YY.LE.O.O))  THEN 
FX  ■  O.O* 

ELSE 

IF (( (ALPHA. GE.O .9) .AND. (ALPHA. LE. 1 . 1)) .OR. 

1  ((BETA. GE. 0.9) .AND. (BETA. LE. 1.1)))  THEN 
FAC1  =  EXP ((ALPHA  -  1 .0)*AL0G(YY) ) 

FAC2  «  EXP ((BETA  -  1 ,0)*AL0G(1.0  -  YY)) 

FX  *  FAC1*FAC2 
ELSE 

FX  =  (YY** (ALPHA  -  1.0))*(1.0  -  YY)**(BETA  -  1.0) 

ENDIF 
END  IF 
RETURN 
END 
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SUBROUTINE  OPTIONS 


C 

C  SUBROUTINE  OPTIONS  ENABLES  THE  INTERACTIVE  USER  OF  BETAFACT  TO 

C  MODIFY  THE  PROBLEM  BEING  EXECUTED  BY  THE  PROGRAM.  THE  SUBROUTINE 

C  LISTS  THE  AVAILABLE  OPTIONS  FOR  THE  USER  AND  PROMPTS  FOR  AN 

C  OPTION  SPECIFICATION.  CONTROL  IS  THEN  RETURNED  TO  THE  MAIN 

C  PROGRAM  WHICH  INTERPRETS  THE  USER  SELECTED  OPTION  AND  PROCEEDS 

C  ACCORDINGLY. 

C 

COMMON  /BLKOl/  IPARAM,  ITYPE ,  ISIZE,  I SEED,  IALGO ,  IXVAL,  ITABL, 

1  IXALG,  IOPT,  NINTV,  NINC,  NINP,  NOUT,  NPLT,  NPRT 

WRITE (NOUT, 1000) 

1000  FORMATC/,’  SELECT  OPTION  TO  qUIT  OR  TO  RUN  A  MODIFIED  VERSION', /, 


1 

t 

OF  THE 

CURRENT  PROBLEM:' ,/,/, 

2 

i 

0 

- 

TERMINATE  PROBLEM.',/, 

3 

) 

1 

- 

CHANGE  INDEPENDENT  PARAMETER  VALUES.',/, 

4 

> 

2 

- 

CHANGE  NUMBER  OF  ALGORITHM  EVALUATIONS',/, 

5 

) 

AND/OR  THE  RANDOM  NUMBER  SEED  INTEGER.',/, 

6 

’ 

3 

- 

CHANGE  PARAMETER  NOMINAL  VALUES  AND/OR',/, 

7 

) 

ASSOCIATED  K-FACTOR  UNCERTAINTIES.',/, 

8 

) 

4 

- 

CHANGE  UNCERTAINTY  DISTRIBUTION  TYPE.',/, 

9 

j 

5 

- 

OPTIONS  1  AND  2.',/, 

* 

) 

6 

- 

OPTIONS  1  AND  3.',/, 

* 

i 

7 

- 

OPTIONS  1  AND  4. ’,/, 

* 

) 

8 

- 

OPTIONS  2  AND  3. ',/, 

* 

> 

9 

- 

OPTIONS  2  AND  4. ',/, 

* 

) 

10 

- 

OPTIONS  3  AND  4. ',/, 

* 

; 

11 

- 

OPTIONS  1,  2,  AND  3. ' ,/, 

* 

) 

12 

- 

OPTIONS  1,  2,  AND  4. ' ,/, 

* 

) 

13 

- 

OPTIONS  1,  3,  AND  4. ’ ,/, 

* 

) 

14 

- 

OPTIONS  2,  3,  AND  4. ' ,/, 

* 

) 

15 

- 

OPTIONS  1,  2,  3,  AND  4. ' ,/, 

* 

) 

16 

- 

COMPLETELY  MODIFY  THE  CURRENT  PROBLEM.’,/) 

READ (NINP, 

*)  IOPT 

RETURN 

END 
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SUBROUTINE  FNRN (NUMNR , RN) 


C 

C  SUBROUTINE  FNRN  RETURNS  A  NORMALLY  DISTRIBUTED  RANDOM  NUMBER 

C  WITH  VARIABLE  NAME  RN.  THE  INTEGER  NUMNR  IS  USED  ONLY  DURING 

C  THE  FIRST  CALL  TO  THE  SUBROUTINE  AND  ITS  VALUE  IS  NOT  CRITICAL. 

C  HOWEVER,  IF  A  DIFFERENT  STRING  OF  RANDOM  NUMBERS  IS  DESIRED 

C  EITHER  DURING  THE  COURSE  OF  A  SINGLE  BETAFACT  ANALYSIS  OR  FROM 

C  ONE  ANALYSIS  TO  ANOTHER,  THE  VALUE  FOR  NUMNR  MUST  BE  CHANGED 

C  AT  THE  BEGINNING  OF  EACH  ANALYSIS.  THE  SUBROUTINE  CALLS 

C  FUNCTION  GASDEV  WHICH  COMPUTES  A  VALUE  FOR  RN. 

C 

DATA  INFIRST  /O/ 

IF  ( INFIRST. EQ.O)  THEN 
INFIRST  *  1 
IF  (NUMNR. GE.O)  THEN 
IDUM  *  -NUMNR 
ELSE 

IDUM  *  NUMNR 
ENDIF 
ENDIF 

RN  *  GASDEV (IDUM) 

RETURN 

END 
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SUBROUTINE  UPR1(NUMUR,RN) 


C 

C  SUBROUTINE  UPR1  RETURNS  A  UNIFORMLY  DISTRIBUTED  RANDOM  NUMBER 

C  WITH  VARIABLE  NAME  RN.  THE  INTEGER  NUMUR  IS  USED  ONLY  DURING 

C  THE  FIRST  CALL  TO  THE  SUBROUTINE  AND  ITS  VALUE  IS  NOT  CRITICAL. 
C  HOWEVER,  IF  A  DIFFERENT  STRING  OF  RANDOM  NUMBERS  IS  DESIRED 
C  EITHER  DURING  THE  COURSE  OF  A  SINGLE  BETAFACT  ANALYSIS  OR  FROM 

C  ONE  ANALYSIS  TO  ANOTHER,  THE  VALUE  FOR  NUMUR  MUST  BE  CHANGED 

C  AT  THE  BEGINNING  OF  EACH  ANALYSIS.  THE  SUBROUTINE  CALLS 

C  FUNCTION  RANI  WHICH  COMPUTES  A  VALUE  FOR  RN. 

C 

DATA  IUFIRST  /O/ 

IF  ( IUFIRST. EQ.O)  THEN 
IUFIRST  =  1 
IF  (NUMUR. GE.O)  THEN 
IDUM  *  -NUMUR 
ELSE 

IDUM  ■  NUMUR 
ENDIF 
END  IF 

RN  =  RANI (IDUM) 

RETURN 

END 
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FUNCTION  GASDEV(IDUM) 


C 

C  THE  FOLLOWING  ROUTINE  IS  TAKEN  FROM: 

C  NUMERICAL  RECIPES  -  THE  ART  OF  SCIENTIFIC  COMPUTING 

C  W.  H.  PRESS,  B.  P.  FLANNERY,  S.  A.  TEUKOLSKY,  AND 

C  W.  T.  VETTERLING;  CAMBRIDGE  UNIVERSITY  PRESS  1988. 

C  P.  202-203. 

C 

C  RETURNS  A  NORMALLY  DISTRIBUTED  DEVIATE  WITH  ZERO  MEAN  AND 

C  UNIT  VARIANCE,  USING  RANl(IDUM)  AS  THE  SOURCE  OF  UNIFORM 

C  DEVIATES . 

C 

DATA  ISET  /O/ 

IF  (ISET.Eq.O)  THEN 
1  VI  =  2 . *RAN1 (IDUM)  -  1. 

V2  =  2 . *RAN1 (IDUM)  -  1. 

RRR  =  Vl**2  +  V2**2 

IF  (RRR.GE.l.)  GO  TO  1 

FAC  =  SQRT(-2.*AL0G'RRR)/RRR) 

GSET  =  V1*FAC 
GASDEV  =  V2*FAC 
ISET  *  1 
ELSE 

GASDEV  =  GSET 
ISET  =  0 
END  IF 
RETURN 
END 
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FUNCTION  RANl(IDUM) 


C 

C  THE  FOLLOWING  ROUTINE  IS  TAKEN  FROM: 

C  NUMERICAL  RECIPES  -  THE  ART  OF  SCIENTIFIC  COMPUTING 

C  W.  H.  PRESS,  B.  P.  FLANNERY,  S.  A.  TEUKOLSKY,  AND 

C  W.  T.  VETTERLING;  CAMBRIDGE  UNIVERSITY  PRESS  1988. 

C  P.  196-197 

C 

C  RETURNS  A  UNIFORM  RANDOM  DEVIATE  BETWEEN  0.0  AND  1.0. 

C  SET  IDUM  TO  ANY  NEGATIVE  VALUE  TO  INITIALIZE  OR 

C  REINITIALIZE  THE  SEQUENCE. 

C 

C  CONSTANTS  ARE  FOR  COMPUTERS  WHICH  OVERFLOW  AT  POSITIVE 

C  INTEGER  VALUES  OF  2**31  -  1. 

C 

DIMENSION  RR(97) 

PARAMETER  (Ml=259200 , IA1=7141 , IC1 =54773 ,RM1=1 . /Ml) 
PARAMETER  (M2- 134456, IA2=8121 ,IC2=28411 ,RM2=1 . /M2) 
PARAMETER  (M3=243000 , IA3=4561 , IC3=51349) 

DATA  IFF  /O / 

IF  (IDUM.LT .0 . OR. IFF . EQ . 0)  THEN 
IFF  =  1 

1X1  =  M0D(IC1-IDUM,M1) 

1X1  =  M0D(IA1*IX1+IC1,M1) 

1X2  =  M0D(IX1 ,M2) 

1X1  =  M0D(IA1*IX1+IC1 ,M1) 

1X3  =  M0D(IX1 ,M3) 

DO  11  J  =  1,  97 

1X1  =  M0D(IA1*IX1+IC1 ,M1) 

1X2  =  M0D(IA2*IX2+IC2 ,M2) 

RR( J)  =  (FLOAT(IXl)  +  FL0AT(IX2)*RM2)*RM1 
11  CONTINUE 
IDUM  =  1 
END  IF 

1X1  =  M0D(IA1*IX1+IC1 ,M1) 

1X2  =  M0D(IA2*IX2+IC2 ,M2) 

1X3  =  M0D(IA3*IX3+IC3 ,M3) 

J  =  1  +  (97*IX3)/M3 
IF  (J.GT.97.QR. J.LT.l)  PAUSE 
RANI  =  RR( J) 

RR( J)  =  (FLOAT(IXl)  +  FL0AT(IX2) *RM2)*RM1 

RETURN 

END 
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